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ABSTRACT 


This  papery  which  is  Based  on  the  technique  published  by  Joel  N. 
Franklin ^gTves  a method  for  generating  digitally  a time  series  with  a 
given  power  spectral  density  function.  A computer  program  to  carry  out 
this  method  was  written  for  the  Control  Data  Corporation  3200  computer, 
and  the  program  was  tried  on  some  specific  example  problems.  Then  for 
each  of  these  examples,  the  power  spectral  density  function  of  the  gen- 
erated time  series  was  compared  with  the  specified,  theoretical  power 
spectral  density  function. 
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I . INTRODUCTION 


A stationary  Gaussian  random  process,  x(t),  is  called  a 

time  series  if  t represents  discrete  values  of  time.  The  collection 

of  finctions  of  time  x(t)  = x^(t)  is  called  a random  process  if  r is 

a random  variable  ranging  over  some  measure  space,  R.  The  process  is 

called  Gaussian  if,  for  every  finite  collection  of  times  t^  < — < t^, 

the  random  variables  x(t^), , x(tn)  have  a multivariate  Gaussian 

distribution.  The  process  is  called  stationary  if,  for  any  increment 

At,  the  random  variables  xCt^+At)  have  the  same  joint  distribution  as  the 

random  variables  xft.).1 

i 

A common  method  for  obtaining  a time  series  is  to  record  on 
magnetic  tape  the  amplitude  of  a physical  quantity,  such  as  the  pressure 
waves  in  a large  body  of  water,  the  atmosphere,  or  the  earth  or  the  waves 
on  the  surface  of  a large  body  of  water.  This  analog  data  can  be  thought 
of  as  a time  series  which  can  be  either  analyzed  by  means  of  an  analog 
computer  or  converted  into  digital  data  and  then  analyzed  by  means  of 
a digital  computer.  This  method  for  obtaining  a time  series  can  be  time 
consuming  and  the  conditions  which  are  necessary  to  obtain  a time  series 
with  specified  statistical  properties  are  often  difficult  to  achieve; 
therefore,  a method  is  desired  for  computing  numerically  a time  series 
with  given  statistical  properties. 

2 

This  paper  summarizes  a technique,  proposed  by  Joel  N.  Franklin  , 
for  computing  a stationary  Gaussian  random  process  with  a given  power 
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spectral  density  function.  It  also  contains  a computer  program  for  the 
Control  Data  3200  computer  which  applies  this  technique  in  computing  a 
specified  stationary  Gaussian  random  process  along  with  four  example 
problems  which  have  been  tried  on  this  program. 

A.  Applications 

The  method  given  in  this  paper  can  be  used  to  generate  a 
stationary  Gaussian  statistical  process  with  a specified  power  spectral 
density  function  which  satisfies  the  conditions  (1.8).  Forms  for  the 
power  spectral  density  function  of  many  statistical  processes  can  be 
found  in  the  literature;  for  example,  the  general  form  for  the  power 
spectral  density  function  of  low  frequency  atmospheric  turbulence  is 
given  by  Bendat.^  Problem  four  in  Chapter  III  gives  a specific  example 
of  such  a power  spectral  density  function. 

As  stated  previously,  a time  series  can  be  obtained  by  recording 
the  amplitude  of  a statistical  process,  such  as  a pressure  wave,  at 
regular  time  intervals.  However,  when  a time  series  is  obtained  in  the 
field,  the  time  series  which  is  recorded  on  magnetic  tape,  called  the 
total  time  series,  is  a combination  of  the  time  series  obtained  from  the 
signal  and  the  time  series  obtained  from  the  noise.  A signal  is  a 
detectable  physical  quantity  or  impulse  by  which  messages  or  information 
can  be  transmitted,  and  noise  is  an  unwanted  detectable  physical  quantity 
or  impulse  which  exists  either  by  itself  or  in  interference  with  a 
specified  signal.  Now  if  the  power  spectral  density  function  of  either 
of  the  three  time  series  mentioned  above  is  known  and  satisfies  the 
conditions  (1.8),  a stationary  Gaussian  random  process,  which  has  the 
same  power  spectral  density  function  as  the  time  series  obtained  from 
the  physical  wave,  can  be  generated  digitally. 
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A problem  of  much  interest  is  the  study  of  detecting  a known 
signal  in  various  noise  backgrounds.  By  the  method  given  in  this  paper 
a time  series  can  be  generated  which  has  the  same  power  spectral  density 
function  as  the  time  series  obtained  from  a specified  noise  background. 

For  example,  suppose  the  time  series  of  a noise-free  signal  of  the 
infrasound  (pressure  waves  having  a frequency  lower  than  about  16  cycles 
per  second)  produced  by  a thunderstorm  is  obtained.  Then  this  time  series 
can  be  viewed  under  any  specified  background  conditions  if  the  power 
spectral  density  function  of  the  time  series  of  the  specified  noise  back- 
ground is  known.  This  can  be  accomplished  by  combining  these  two  time 
series  with  the  use  of  either  the  digital  or  analog  facilities  of  a 
computer.  The  magnitude  of  the  time  series  of  the  signal  can  be  varied 
in  relation  to  the  magnitude  of  the  time  series  of  the  noise,  and  by  this 
method,  the  time  series  of  the  infrasound  produced  by  a thunderstorm  at 
a specified  distance  with  a specified  background  can  be  simulated.  The 
power  spectral  density  functions  of  these  total  time  series  can  be  used 
in  developing  methods  for  detecting  and  tracking  thunderstorms  by  giving 
examples  of  the  power  spectral  densities  of  the  infrasound  produced  by 
thunderstorms  at  various  distances  with  various  backgrounds . 

The  ability  to  obtain  a time  series  for  a given  power  spectral 
density  function  could  also  be  useful  in  the  advancement  of  the  theory 
of  optimum  linear  prediction  and  filtering.  This  theory  is  used  in 
designing  engineering  systems  which  either  project  into  the  future  by 
information  obtained  in  the  past  or  recover  desired  signals  which  have 
been  distorted  by  random  noise  disturbances.  These  systems  can  be  applied 
to  communication,  meteorological  forecasting,  and  economic  analysis. 

A specific  example  is  the  technique  of  combining  two  independently 
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obtained  sets  of  related  perturbed  messages  to  obtain  an  improved  estimate 

4 

of  the  message. 

The  last  application  to  be  mentioned  is  the  study  of  smoothing 

5 

techniques  such  as  the  harming  and  hamming  windows.  Here  it  is  of  interest 
to  learn  how  the  minor  lobes  of  a power  spectral  density  function  are 
affected  by  the  major  lobe  when  these  smoothing  techniques  a-'e  applied. 

This  can  be  done  by  constructing  time  series  for  which  the  distances  be- 
tween the  major  and  minor  lobes  of  their  power  spectral  density  functions 
vary.  The  relative  size  of  these  lobes  can  also  be  varied  along  with  the 
distance  between  them;  therefore  the  effect  of  a smoothirg  technique  on  a 
given  power  spectral  density  function  can  be  studied. 

B.  Summary  of  Method 

The  method  under  consideration  for  computing  a stationary 
Gaussian  random  process  is  the  method  proposed  by  Joel  N.  Franklin. ^ 

In  order  to  compute  such  a random  process  by  this  method,  either  the 
autocorrelation  function,  which  is  defined  in  terms  of  the  expected  value 
operator,  or  the  power  spectral  density  function  must  be  known. 

The  expected  value  operator  E is  defined  as 


IN 

4i]  = X: 


d.i) 


if  the  limit  exists.  For  a finite  sequence  x , the  operation  lim  is 

1 N-»  w 

neglected.  Throughout  this  paper  it  is  assumed  that  E[x(t^)]=0  for  each 
time  series  x(t). 


The  autocorrelation  function  is  defined  as 


R(ti>t2)  = E[x(tL)x(t2)] 


which  is  written  as 


(1.2) 


R(t)  a E[x(t)x(t-T)] 


(1.3) 


t 
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if  x(t)  is  a stationary  process. 

The  power  spectral  density  function  for  a stationary  process  is 

defined  as 


if  the  integral  exists. 

The  Fourier  transform  indicated  by  (l.4)  can  generally  be  inverted; 
therefor*  , if  V(ui)  is  known  and  can  be  inverted,  R(t)  can  be  found.  Since 
it  is  more  natural  to  work  with  the  power  spectral  density  function  and 
since  Vfto)  and  R(t)  are  generally  interchangeable  this  method  is  designed 
to  construct  a stationary  Gaussian  random  process  from  a given  power  spec- 
tral density  function.  Even  though  V(co)  and  R(t)  are  interchangeable,  small 
perturbations  in  V(<jo)  may  cause  great  changes  in  R(t).  This  case  is  illus- 
trated by  example  problem  four  in  Chapter  III. 

Franklin  attacks  the  problem  by  considering  the  linear  trans- 
formation 

x(t)  =»  j'  g(t-s)v(s)ds  . (1.5) 

— 00 

He  then  sets  V equal  to  the  Dirac  delta  function;  therefore  g(t)  is  the 
response  of  a filter  to  a delta  function  input.  It  is  shown  by  Davenport 

•J 

and  Root  that  the  power  spectral  density  functions  vx(“>)  and  Vv(cu)  of 
x(t)  and  v(t),  respectively,  are  related  by 

Vx(u>)  = | G(u>)  | 2 Vv  (u>)  , (1.6) 

where  G(ou)  is  the  Fourier  transform  of  g(t). 

Now  if  a random  process  v(t),  known  as  white  noise  which  has  the 
property  that  Vv(u>)=l,  is  obtained,  (1.6)  becomes 


I 


V(o>)  = Vx(o>)  = | G(a>)  | 

g 

It  is  shown  by  Davenport  and  Root  that  if 
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(1.7) 


0<V(co)<°°,  V(cjd)  = V(-ca),  and  V(u>)-*0  as  or*±°°,  (1.8) 

and  that  if  V(o))  can  be  represented  as  a rational  function  (i.e.,  the 
quotient  of  two  polynomials)  in  oo,  then  for  real  u>,  V(cjd)  can  be  represented 
by 

where  P and  Q,  are  polynomials  with  real  coefficients  and  of  degree  m and  n, 
respectively,  m<n,  and  the  zeros  of  Q(|)  lie  in  the  half-plane  Re  £<0.  Now 
by  (1.7)  and  (1.9),  an  obvious  choice  for  G(oo)  is 

aM  - IBs}  - (i-io) 

and  if  D is  the  differential  operator,  d/dt,  we  have 

x(t)  = Iff}  v(t)  * (l-ll) 

The  specified  time  series  x(t)  is  found  by  first  finding  a steady-state 
solution  to  the  differential  equation 

Q(D)cp(t)  = v(t)  , (1.12) 

for  all  real  t and  then  combining  the  derivatives  of  <p(t)  of  order  less 
than  n by 

x(t)  = P(D)qKt)  . (1.13) 

The  first  step  in  finding  a solution  <p(t)  to  the  differential 
equation  (1.12)  is  to  generate  a completely  equidistributed  sequence 
which  is  also  a "white  sequence."  This  sequence  is  to  be  used  as  the 
digital  representation  of  a white  noise  source. 

A "white  sequence"  is  defined  to  be  a sequence  for  which 

R(t)  = 0 for  t * 0 . (1.1*0 

A sequence  is  said  to  be  equidistributed  by  k's  if  k is  a 
positive  integer  and  if  for  every  set  of  k intervals  (ai/b^),  i=l,  , k. 
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4 
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where  0 s a^<  K S 1, 


*1 


(b^a.^)  as  N -»  <*>  , 


(1.15) 


a,  £ v , .<  b.  i=l 
i n+i  i 

(i-1,— ,k) 

n*l,  — ,N 

and  a sequence  equidistributed  by  k's  for  all  k is  called  a completely 
equidistributed  sequence. 

Let  the  notation  (cp)  denote  the  fractional  part  of  Cp.  It  is 
shown  by  Franklin  that  for  almost  every  0 > 1,  (011}  is  a completely  equi- 
distributed sequence^1;  furthermore  he  suggests  that  the  sequence  Vn=(0n), 
where  9 is  a transcendental  number  greater  than  one , be  used  as  the  required 
"white  sequence."  For  example,  if  6 = 2.72,  » (2.72)  = .72, 

v2  * {2.722}  = .3984,  etc. 

Because  of  computational  difficulties  which  limit  the  length  of 


(1) 


(s) 


for  a given  9 > 1,  several  different  sequences,  , , , 

(n  = 0,1,2, ),  may  be  needed  for  a single  application.  These  sequences 


should  be  interlaced  by 

_ M 


V = V' 

ns+1  n 


’ Vns+2 


v(2)  ... 

Vn  ' 


v = v''5',  n = 0,1,2,—  (1.16) 
ns+s  n ’ ’ 


= v(s) 
n 

to  form  the  required  "white  sequence."  A few  of  these  difficulties  will  be 

pointed  out  in  Chapter  II,  Section  A. 

Now  that  an  equidistributed  sequence  has  been  constructed  on 

the  range  0 § v < 1,  a sequence  w^,  which  will  be  called  the  w-sequence,  of 

independent  samples  from  the  Gaussian  distribution  with  mean  0 and  vari- 
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ance  1 can  be  constructed  by  the  method  of  Box  and  Muller, 
is  given  by  the  two  formulas, 

W2n-1  = ("2  10  v2n-l}^  COS  2nVS 


This  method 


r2n 


2n 


= (-2  In  v„  sin  2nv„  for  n=»l,2,3,— 


(1.17) 


2n-l 


2n 


8 


L 
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Now  that  the  "white  sequence,"  v^,  and  the  w-cequence,  w , have 
"been  defined,  the  next  step  in  finding  a solution  to  the  differential 
equation  (1.12)  is  to  consider  the  given  power  spectral  density  function 
V(u>),  which  must  he  a rational  function  satisfying  the  conditions  (1.8), 
and  the  given  time  increment  At  at  which  samples  of  the  specified  time 
series  are  to  he  calculated. 

The  equation  (1.12)  can  he  represented  as 

v(n)(t)  + «2.p(n-1>(t) 


+ — - + an+^<p(t)  = v(t)  for  -«<t<oo  (l.l8) 


where  a^=l  and  a.,  ISiSn+l,  are  the  real  coefficients  of  the  polynomial  Q. 


The  vectors 


:(t)  = 


zL(t) 

z2(t) 


Zn(t) 


<P(t) 

cp'(t) 


, t = 0, At,  2At, 


(1.19) 


must  he  found  in  order  to  compute  the  specified  time  series  hy  the  equation 
(1.13).  The  vector  z satisfies  the  differential  equation 

= A z(t)  + f(t)  , -CO  < t < 00  (1.20) 

where 

r 

0 1 0 • • 0 

0 0 1 • • 0 

1 

(1.21) 


0 0 0 • • 1 

■fit  , . 

n+l  n n-1  2 


, ..dm* 


(1.22) 


In  order  to  find  z(0),  the  solution  vector  (m  m , m ) must 

" V 2 n 


first  be  found  to  the  set  of  n linear  equations  in  n unknowns  given  by 


£<■ 


t + is±|iihi 


.q-1  _ |0,  k=l,  n-1 

an-2(q-l)+kmq  “u/2,  k=n 


(1.23) 


The  elements  of  this  solution  vector  are  then  placed  into  the  matrix  M by 

0,  i+j  odd 

mii  = \ (-l)^-i^/2n,  s ( 1 -2' 

J ' m(j+i)/2  , i+j  even 


(1.24) 


for  1 § i,  j § n. 


The  formulas  (1.23)  and  (1.24)  are  proved  by  Franklin. ^ M is  the  positive 

definite  moment  matrix  for  the  vector  s(t)  for  all  time  t,  including  t = 0. 

Since  M is  a positive  definite,  real,  symmetric  matrix,  there  exists  a 

lower-triangular  matrix  T with  positive  elements  on  the  main  diagonal  such 
* 12 

that  M = T T. 


Now  it  becomes  necessary  to  define  the  sequence  of  n-dimensional 


vectors, 


(1.25) 


i 
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whose  elements  are  the  elements  of  the  w-sequence,  and  z(0)  is  determined 


z(0)  =>  T w^U;  . (1.2 1 

Franklin  shows  that  if  z(kAt)  is  known,  the  following  steps  are 

needed  to  compute  z((k+l)At).  z(0)  has  been  computed  in  equation  (1.26); 

therefore,  it  is  used  as  the  starting  value  and  z((k+l)At)  is  found  by 

induction.  The  first  step  in  finding  z((k+l)At)  is  to  compute  the  matrix 

exp(AtA)  with  elements  e. .,  1 S i,  j § n. 

i J 

It  then  becomes  necessary  to  solve  the  set  of  simultaneous 


equations 


lb.  ,=e . e.  , i<n  or  j<n 
\ I ij  in  jn 

u + a u ) =</'  ° „ u 

lk^kj  jk^k;  V 2 . 

0 ° |b  =e  -1  i=n  and  j-n 

nn  nn 


where  u. . are  the  unknowns  and  a. . are  the  elements  of  A for  1 5 i,j  S n. 

1.1  ij 

The  elements  p are  also  the  elements  of  the  symmetric  moment  matrix  M ; 

^ ~ r 

there fore, u.  . = u,.  . It  also  can  be  seen  that  b.  . = b.. ; therefore,  the 
' ij  Ji  iJ  Ji 

2 2 

n equations  in  n unknowns  represented  by  (l.27)  can  be  reduced  to  the 


^■n(n+l)  equations  in  £n(n+l)  unknovms, 


x*, 


ik^jk  + 


aik^kj 


a0k^ik  + , 


\jk^ki  ^ij 


for  1 S j S i S n. 

When  the  solution  to  (1.28)  has  been  found,  the  positive  definite,  real, 

symmetric  matrix  Mr  is  known,  and  therefore  a lower  triangular  matrix  T^ 

* 

can  be  found  such  that  M = T T 

r r r 

(k) 

T^,  At,  exp(AtA),  z(0),  and  wv  , k=l,2, , are  now  known,  and 

the  last  step  in  evaluating  z((k+l)At)  is  carried  out  by 

z((k+l)At)  = exp(AtA)z(kAt)  + T w^+1)  (l.2S 

r 

for  k = 0, 1, 2, — 


11 

Now  the  solution  <p(t)  to  the  equation  (1.12)  and  its  first  (n-l) 
derivatives  are  known.  Equation  (I.13)  becomes 

x(t)  = 'b1zm+1(t)  + bgzm(t)  + — + bm+1zL(t),  t = 0, At, 2At, (1.30) 

where  b^  # 0 and  b^,  1 s i S m+1,  are  the  real  coefficients  of  the  poly- 
nomial P and  z^(t),  1 S i S n,  are  the  elements  of  the  vector  z(t)  as 
defined  by  (1.19).  This  completes  the  process  since  the  specified  time 
series  x(t)  has  been  obtained. 


I 

I 
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II.  COMPUTER  PROGRAMS 


The  method  for  calculating  a time  series  with  specified  power 
spectral  density  function,  as  presented  in  Chapter  I,  is  divided  into 
three  separate  programs.  The  first  program,  ISLASHO  with  subroutine 
RANDOM,  generates  a "white  sequence"  for  a given  d and  buffers  this 
sequence  out  onto  tape.  The  program  CLLCNVRT,  pronounced  Call  Convert, 
with  subroutine  CONVERT  then  converts  the  "white  sequence"  into  the 
w-sequence  by  (1.17).  This  w-sequence  is  then  buffered  out  onto  a 
second  tape  which  is  the  input  tape  to  the  program  GAUSSIAN.  GAUSSIAN 
then  shapes  the  w-sequence  into  a time  series  with  a given  power  spectral 
density  function.  The  power  spectral  density  function  is  defined  in  this 
program  by  the  coefficients  of  the  polynomials  P and  Q as  defined  by  (1.9). 
This  specified  time  series  is  buffered  out  onto  a third  tape.  The  same 
w-sequence  can  be  shaped  into  an  specified  time  series. 

A.  Generating  a "White  Sequence" 

The  purpose  of  the  program  ISIASHO  is  to  generate  a "white 
sequence"  as  described  in  Chapter  I,  Section  B. 

For  the  purpose  of  this  paper,  the  transcendental  number 
greater  than  one  was  chosen  to  be  e = 2. 71828103.  The  "white  sequence", 
v,  = ten)  for  1 S n S 4266,  was  then  generated  on  the  Control  Data  320 0 
computer. 

The  following  approach  was  taken  to  the  problem  of  generating 
a "white  sequence."  The  approach  It.  -ds  to  a description  of  the  method 
used. 

For  every  n > 1,  en  can  no  longer  be  stored  to  full  accuracy  as 
as  an  eleven  digit,  floating  point  number;  therefore,  this  is  obviously 
not  the  approach  to  take. 
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It  can  be  seen  that  about  17^0  digits  are  necessary  to  store  the 
Uooo  4q00 

integer  part  of  e since  log  (e  ) ® 1740.  Initially  e is  rounded  to 
nine  digits;  then  for  the  assurance  that  the  fractional  part  of  en,  v^, 

will  not  start  repeating,  en  is  saved  with  no  round-off  error.  Since  the 

present  approximation  of  e contains  eight  digits  to  the  right  of  the  decimal 

1*000  1^000 
point,  (e  } contains  32,000  digits,  and  therefore  e consists  of 

approximately  33,7^0  digits.  Now  if  these  digits  were  stored  in  groups  of 

4ooo 

threes  as  fixed  point  numbers,  e would  require  about  12,000  fixed  point 
numbers  for  storage.  The  storage  capacity  of  the  Control  Data  3200  is  about 
32,000  24-bit  words;  therefore,  if  two  12,000-word  arrays  are  needed  for  the 
multiplication,  the  machine  still  has  tie  capacity  to  handle  the  program. 

The  largest  fixed  point  number  that  can  be  stored  in  a 24- bit  word 
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is  2 -1  * 8,388,000;  therefore,  all  six-digit  and  most  seven-digit  integers 

can  be  stored  in  one  24-bit  word. 

Now  the  properties  that  the  Control  Data  3200  can  handle  two 
12,000,  24-bit  word  arrays  and  all  six-digit  integers  can  be  stored  in  a 
24-bit  word  will  be  used  in  the  program  ISLASHO  to  raise  any  seven-,  eight-, 
or  nine-digit  number  to  powers  and  save  the  fractional  part.  The  two 
12,000-word  arrays,  IX  and  IY,  are  used  to  store  the  products  as  three-digit, 
fixed  point  numbers.  These  products  should  be  thought  of  as  numbers  repre- 
sented in  the  base  1000.  The  notation  used  in  equations  (2.l),  (2.4),  and 
(2.7)  is  that  of  "long  hand"  multiplication. 

The  number  which  is  to  be  raised  to  powers  is  read  into  ISLASHO 
as  the  three,  three-digit,  fixed  point  numbers,  1X3,  1X2,  and  1X1,  with  1X1 
containing  the  three  right-most  digits,  an  the  method  used  in  ISLASHO,  we 
first  set  IX(  1 ) = 1X1,  IX(2)  «•  1X2,  IX(3)  = 1X3,  and  L = 3 


and  then  start  the  following  procedure. 


In  step  one  we  have 
IX(L)  •••• 
IY(L+l)  IY(L)  


IX(3)  IX(2)  IX(1) 

1X1 

TyT31  tyT21  iyTT) 


(2.1) 


IY(l)  is  computed  by 


(2.2) 


also  let 


ICABRY1  = I 


IX(l)  X 1X1 


1000 


R j~|  is  defined  to  be  the  remainder  of  a/b,  and  I |-gj  is  defined 

to  be  the  integral  number  of  times  b divides  a where  a and  b are  integers. 

Then  for  2 5 i S L, 

let  IY( i ) = R(T. ) and  ICARRY.  = l(T. ) 

l i i 

where 

flX(i)  x 1X1  + ICARRY.  _A 


T. 

l 


1000 


(2.3) 


To  terminate  this  step,  set  IY(L+l)  = ICARRY^  and  increase  L by  one. 

In  step  two  we  have 

IX(L-l)  •••  IX(3)  IX( 2)  IX(1) 

IX  2 

IY(L+l)  TYTlI  IY(L-1)  •••  IY( 3)  IY(2)  IY( l)  * 


(2.4) 


where  IY(l)  has  the  same  value  as  in  step  one,  L-l  has  the  same  integer 


value  as  L in  step  one,  and  TY(2)  = R(T^)  ICARRY^  = I(T^)  for 

T _ IX(1)  X 1X2  + IY( 2 ) 
r2  = 1000 

Let  IY( i)  denote  IY(i)  as  computed  in  the  previous  step. 


(2.5) 


where 


Let  IY( i ) = R(T.)  and  ICARRY ^ = l(Ti), 


IX(i-l)  X 1X2  + IY( i ) + ICARRY._l 
1000 


for  ) i i i L 


(2.6) 
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r 


i 

i 

i 


Step  two  is  also  terminated  by  setting  IY(L+l)  = ICARRY^  and  then  increasing 
L by  one. 


In  step  three  we  have 

IX( L-2) 


IX(3)  IX( 2)  IX(1) 

2L 


IY(L+1)  IY[l1  IY(L-l)  IY(L-2)  •••  iY(3)  u{2)  iy(I) 

where  IY(l)  and  IY(2)  remain  invariant  from  step  two,  IY(3) 

ICARRY^  = I(T^)  where 

_ / ix( 1 ) x 1X3  + iy(1T 
3 " iooo 


(2.7) 


R(Tj),  and 


Let  IY( i)  = R(T. ) and  ICARRY.  = I (T. ) 


where 


T.  = 
i 


for  4 5 i s l. 


IX(i-2)  x 1X3  + IY(i)  + ICARRY. 


i-1 


1000 


(2.8) 


(2.9) 


To  terminate  step  three  we  set 

|iCARRYL,  if  ICARRYl  * 0 
IY(L+I)  =<  not  defined,  if  ICARRYl  = 0 

~ _ J L+l,  if  ICARRY  * 0 
L ^L,  if  ICARRYl  = 0 , 


| 


(2.10) 

(2.11) 


and  then  L = L. 

It  can  be  seen  from  equations  (2.l)  through  (2.y)  that  no  fixed 

2 

point  number  calculated  by  the  above  method  will  be  greater  than  999  + 999+ 

999  = 999>999  which  is  a six-digit  number  and  can  therefore  be  bandied  by 


the  Control  Data  3200. 

The  IY  array  is  then  placed  into  the  IX  array  said  eight,  nine, 
or  ten  correct  digits  (depending  on  where  the  theoretical  decimal  point 
is  located  in  relation  to  the  three  digits  of  IX(i)  which  contains  the 
theoretical  decimal  point)  after  the  theoretical  decimal  point  are  placed 
into  the  floating  point  X array.  The  X array  therefore  contains  numbers 
between  zero  and  one;  this  is  the  desired  "white  sequence."  The  process, 
starting  with  step  one  and  continuing  through  the  placing  of  a new  number 
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i 

l 

l 

[ 

i 

! 

I 

t 

1 


into  the  X array,  is  repeated  until  L is  even  and.  ISTOP  < L S 12,000  where 
ISTOP  is  an  input  parame’er. 

The  entire  IX  array  is  placed  on  tape  tefore  ISLASHO  terminates; 
therefore,  this  calculation  can  later  he  continued.  However,  in  order  to 
calculate  many  more  than  4266  points  in  this  "white  sequence"  hy  this 
method,  ISLASHO  would  have  to  be  modified  such  that  part  of  the  IX  and  IY 
arrays  could  be  stored  elsewhere  than  in  the  memory  of  the  computer  during 
the  calculation  of  each  point  of  X because  these  two  arrays  would  soon  out 
grow  the  storage  capacity  of  the  computer.  Also  the  longer  these  arrays 
become,  the  time  needed  to  calculate  each  point  of  X increases  until  this 
method  would  no  longer  be  practical.  For  example,  v , 1 S n S 500,  was 
computed  in  about  45  seconds  as  compared  with  thirteen  minutes  needed  to 
compute  the  500  points,  v , 3501  s n § 4000.  It  was  noted  that  if  t minutes 
were  reeded  to  compute  a block  of  500  points  of  v^,  approximately  (t  + 1.75) 
minutes  were  needed  to  compute  the  next  500  points  of  v . It  therefore  may 
become  necessary  to  generate  several  "white  sequences,"  vj^,  •••  V^S^,  for 

different  transcendental  numbers  and  then  interlace  these  sequences  by  the 
method  (l.l6)  to  form  the  desired  "white  sequence." 


I 

I 

I 


I 


I 

I 

I 

[ 


I 

I 


[ 

l 

1 

I 


1 

I 

l 
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The  following  is  the  list  of  input  parameters  to  the  program 
ISLASHO.  The  numbers  in  parentheses  were  the  numbers  used  to  find 
(en)  for  1 5 n § 4266. 


NO  DEC 


1X3,  1X2,  1X1 


IF  PRNT  SQ,  = 1, 


- 0, 


IF  PRNT  ME  = 1, 


ISTOP 


- o. 


IBUF 


the  number  of  digits  to  the  right  of  the  decimal  point 
of  the  number  which  is  to  be  raised  to  power.  (8). 
these  three  parameters  each  contain  three  digits  of  the 
number  which  is  to  be  raised  to  pow-  ’~s.  (2.71828l83  was 
used;  therefore,  1X3  = 271,  1X2  = 828,  and  1X1  = 183). 
the  X array  is  buffered  out  onto  tape  and  printed  in 
blocks  of  IBUF  points, 

the  X array  is  only  buffered  out  on  tape  in  biocks 
of  IBUF  points,  (l). 

buffer  out  onto  tape  and  print  the  entire  IX  array  just 
before  the  program  terminates, 

only  buffer  out  onto  tape  the  entire  IX  array  (0). 
the  program  terminates  when  L is  even  and  ISTOP  < L § 
12,000.  ISTOP  should  never  exceed  six  less  than  the 
dimension  on  the  IX  and  IY  arrays  (11 99° ) - 
the  number  of  points  of  the  "white  sequence,"  X array, 
that  should  be  handled  in  one  block.  IBUF  should  be  sin 
even  number  and  not  exceed  the  dimension  of  X which  is 
500  in  the  listing.  Also  IBUF  £ 388  if  the  subroutine 
SOLUTION  in  GAUSSIAN  is  not  to  be  changed. 


. 

I 


■ 


d 


The  output  tap-  has  the  following  form: 


ICOUNT 


I DEC 


L 


IX(1) 


IBUF 

x(i)  • • • • X(IBUF) 


IBUF 

X(l)  • • • • x(ibuf) 

ICOUNT 

X(L)  . . . . X( ICOUNT) 
end-of-file 

IDEC 

L 

IX(1)  • • • • IX( L) 
end-of-file 

the  number  of  points  in  the  last  block  of  the  X array. 
ICOUNT  S IBUF. 

the  number  of  digits  in  the  decimal  part  of  the  last 
number  in  memory. 

the  number  of  three-digit  numbers  which  compose  the  last 
number  in  memory. 

IX( L)  the  last  number  in  memory  put  on  tape  as  three-digit 
numbers.  IX(l)  is  the  right  most  and  IX(L)  is  the 
the  left  most  three  digits  of  this  number. 


I 
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I 

I 


I 

I 


L 

1 

I 

f 

I 


I 


B.  Constructing  the  w-sequence 

The  only  purpose  of  CLLCNVRT  is  to  convert  a "white  sequence" 
into  its  w-sequence  by  (1.17)  and.  place  this  w-sequence  onto  tape.  This 
program  is  almost  trivial,  hut  it  simplifies  the  programming  of  the  problem. 
It  can  be  seen  by  (1.17)  why  it  is  desired  that  an  even  number  of  points 
of  the  "white  sequence"  be  placed  on  tape  in  each  block.  If  the  number  of 
points  in  a block  is  odd,  the  last  point  in  that  block  will  not  be  used  by 
the  subroutine  CONVERT  in  forming  the  w-sequence. 

The  program  CLLCNVRT  has  one  input  parameter,  IF  PRINT. 

IF  PRINT  = 1,  print  and  buffer  out  the  w-sequence  onto  tape, 

= 0,  only  buffer  out  the  w-sequence  onto  tape. 

The  output  tape  has  the  following  form: 

I BUF 

W(l)  W(I  BUF) 


I BUF 

W(l)  •••  W(I  BUF) 
end-of-file . 

I BUF  same  as  on  input  tape. 

W(l)  •••  W( I BUF)  a block  of  w-sequence  points. 
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C . Constructing  a Specified  Time  Series 

The  program  GAUSSIAN  shapes  the  input  w-sequence  into  a time 
series  x(t)  with  a given  power  spectral  density  function  and  a given 
sampling  time  interval.  The  power  spectral  density  function  must  be 
defined  by  the  real  numbers  a^,  1 S i S n + 1,  and  b^,  1 S i S m+1  for 


P(ico) 

2 

b^(ko)m+  * 

• + bm(ia))  + bm+l 

Q(iu>) 

a-L(iu>)n  + • 

• + an(icw)  + Vl 

v(u>)  = 


where  m < n,  u>  is  real,  and  the  zeros  of  a^i  + 


(2.12) 


+ a I + a . lie  in 
n n+1 


the  halfplane  Re  (?)  < 0.  It  should  be  noted  that  each  of  the  n-dimensional 
vectors  in  (1.25)  is  used  to  construct  one  point  of  the  specified 
time  series;  therefore,  if  the  w-sequence  consists  of  J points,  the 
specified  time  series  will  contain  J/n  points. 

The  numbers  a.,  1 5 i § n+1,  are  read  into  the  vector  SA;  the 

l 

numbers  b^,  1 £ i § m+1,  are  read  into  the  vector  SB;  and  the  time  interval 
At  is  read  into  DT. 

Subroutine  FIND  B places  the  numbers,  a^,  1 £ i § n+1,  into  a 
constant,  n x n matrix  B and  constructs  the  n-dimensional  constant  vector 
[0,  •••,  0,  £]  by  (1.23). 

The  subroutine  GAUSS  P then  finds  the  solution  vector  [m  , • • • m ] 

1 n 

by  means  of  Gaussian  elimination  with  partial  pivoting. ^ Since  the  deter- 
minant of  the  input  matrix  can  be  found  by  simply  computing  the  product  of 
the  elements  on  the  main  diagonal  and  multiplying  this  product  by  the  proper 
sign  after  the  input  matrix  has  been  reduced  to  upper  triangular  form,  the 
determinant  of  the  input  matrix  is  always  calculated  by  GAUSS  P and  printed 
out  by  the  program  GAUSSIAN  as  a guide  to  the  reliability  of  the  solution 


vector. 
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The  subroutine  FIND  CM  places  nu,  1 S i § n,  into  the  positive 
definite  moment  matrix  M whose  elements,  m.  1 s i,  j S n,  have  been 
found  by  (l.2k).  Then  as  stated  in  Chapter  I,  Section  B,  a lower  trian- 
gular matrix  T must  be  found  such  that 

M = T T. * (2.13 

Since  M is  a positive  definite,  real,  symmetric  matrix,  the 
desired  lower  triangular  matrix  T can  be  found  by  Crout  factorization. 

Let  t.  . be  the  elements  of  T and  m.  . be  the  elements  of  M.  It  is  known 
ij 

that  t . = 0 for  j > i ; therefore,  only  t. . for  j § i need  to  be  computed. 

O 1 J 

Then  by  (2.13) 


n.  . = / t. , t .. 
ij  Ui  ik  jk 

k=l 


can  be  obtained.  It  is  then  suggested  that  t. .,  for  j 5 i,  be  computed  in 

— J 

the  following  order: 

ij  = 11,  21,  •••>  nl;  22,  32,  •*,  n2;  •••;  nn. 

By  this  order. 


tll  = mll 


(2.15) 


is  computed  first,  and  the  other  elements  in  the  first  column  are  given 


t. ..  = t. , m. ..  for  i = 2, 
ll  11  ll  ' 


(2.16) 


If  the  preceding  columns  k < j have  been  computed,  the  j — - diagonal  element 


can  be  computed  by 


tij  = mi j 


■1  ‘ttl* 


(2.17) 
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Now  the  elements  below  the  diagonal  can  be  computed  by 

j-1 


'U 


lU 


X 

k=l 


tik  tjk 


(2.18) 


for  i = j+1,  •••,  n where  j < n.  ~ This  method  of  Crout  factorization  is 
carried  out  in  the  subroutine  T T STAR. 

The  matrix  exp(AtA),  where  A has  the  form  as  shown  in  (l.2l)  and 
is  constructed  by  the  subroutine  MAKE  A,  is  computed  by  the  subroutine 
EXP  ADT.  This  exponential  matrix  is  defined  by 


/-, .X  At  . ^ At2  .2  Atk  .k  ^ 

exp(AtA)  = I + -yr  A + — A + • • + y-p-  A + 

which  is  the  same  as 

exp(AtA)  = I + A [ij  + A [yy  a] 

/a*’ 

+ 


(2.19) 


( At\ 

.fat2 

Zl  (&t\ 

A *3 

\ y 

aL— # 

J vv 

A —rr  A 

(2.20) 


...  + ('*’) 
\k  / 


rAtk_1  , k- 

Tk^T:  A 


th 

The  method  of  summation  given  by  (2.20),  where  the  (k-l) — term  is  multi- 
plied by  (At/k)  A and  then  added  to  the  previous  sum,  is  used  to  calculate 
exp(AtA)  in  EXP  A DT. 


The  summation  terminates  when  the  condition  that  b.  .<  c.  max 
, 10  chk 


(Ianil  X 10 


-12 


1 S i § n and  is  an  element  of  A)  holds  for  every 

Atk  k 

b.  .,  1 3 i,  j § n,  where  b.  . are  the  elements  of  B = — — A for  k £ 1.  The 
1 J lj  K 

maximum  of  the  above  set  was  taken  instead  of  the  minimum  so  that  the  computer 
would  not  set  € equal  to  some  very  small  floating  point  numbers  which 
should  actually  be  zero.  However,  if  all  a^  are  approximately  the  same 
in  magnitude,  no  problem  should  arise.  If  vary  greatly  in  magnitude, 
it  may  become  necessary  to  define  a floating  point  zero  and  change  the  sub- 


routine to  check  for  the  minimum  value  of  (la  . I 1 & i S n) . 

ni 
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Let  a,  te  the  elements  of  A where 

id 


A = AxB  = A (~  Ak)  = (~)  AAk  = A*)  A = BxA  . (2.21) 

Now  since  AxB  = BxA,  the  matrix  multiplication  for  (2.2l)  can  be  defined  by 


and 


a.  , = b.  UiS  n-1,  1 S j S n 

id  in,  j' 

A 

a . = b a . 
nl  nn  nl, 

a j * b . + b a 1 S j S n-1 

n, d+1  nd  nn  n,d+l 


(2.22) 


This  is  the  method  of  matrix  multiplication  used  in  EXPADT. 

Now  subroutine  BIGSET  constructs  the  ^ n(n+l)  equations  in  the 

same  number  of  unknowns  by  (1.28)  where  the  b.  . for  1 § j S i § n are 

defined  by  (l.27).  The  unknowns  n , 1 § j 5 i 4 n,  can  be  written  in  the 

1 J 

form  of  an  nxn  lower  triangular  matrix  L,  and  a one-to-one  mapping  f can 
be  defined  which  maps  the  elements  on  and  below  the  main  diagonal  of  L onto 

an  n-dimensional  vector  v where  n = £ n(n+l).  This  vector  is  then  the 

A A 

vector  of  unknowns  associated  with  the  nxn  constant  matrix  defined  by 
(1.28).  After  the  value  of  v has  been  found  by  GAUSS  P,  an  inverse  mapping 
can  be  defined  to  take  the  elements  of  v back  into  L.  L can  then  be  made 
irto  the  desired  real  symmetric  matrix  by  replacing  the  zeros  above  the 
main  diagonal  with  the  corresponding  elements  below  the  diagonal. 

Let  f:L  -*  v be  the  mapping 

0 • • • 0 

11 

^21  ^22  * 


f : 


nl  n2  nn^ 

(nLL>  »2l’  %1;  422'  "**'  4n2’  ^33 

= (Vv2>  •**>  v-)  . 


’ ^nn^ 


(2.23) 


When  ^ is  under  consideration,  all  n elements  in  the  first  j-1 

t-h 

columns  plus  the  upper  i-1  elements  in  the  j—  column  less  the  upper  trian- 
gular set  of  zero  elements  have  been  mapped  into  v;  therefore,  it  can  intui- 
tively be  seen  that  for  1 £ j § i S n,  and  v^,  1 § k 5 n,  the  mapping 

f can  be  defined  by 


which  reduces  to 


k * i + n (j-1)  - [j j (j+l)] 


k = i + (j-l)(2n-j)/2  . 


(2.24) 


(2.25) 


It  must  be  shown  that  f,  as  defined  by  (2.24),  is  the  mapping 
indicated  by  (2.2 5).  Let  f(u. .)  = (v  ) be  denoted  by  f (i,j)  = k for 

1 J K 

1 § j § i § n and  1 £ k § n.  Therefore  show  that  if  1 5 j = j 5 n - 1 and 
i - i = 1,  or  if  i = n,  1 S J s n-1,  and  T = 3 = j+1,  then  k-k  = f(T,3)-f( i, j ) 


Case  I: 


Case  II: 


Let  1 S j = j & n-1  and  i-i  = 1;  this  implies  that 
k-k  = f(I,j)  - f(i, j) 

= I + (j-l)(2n-j)/2-i-(j-l)(2n-j)/2 
= I - i = 1 . 


Let  i=n,  1 5 j 5 n-1,  and  i=j  = j+1;  this  implies  that 
k-k  = f(T,  j)  - f(i,j) 

= I + (j-l)(2n-j)/2  - i -( j-l)( 2n- j )/2 
= j+1  + (j+l-l)(2n-j-l)/2 
-n  - (j-l)(2n-j)/2 
= j+1  + (2nj-j2-j)/2 
-n  - (2nj-j2-2n+j)/2 
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Therefore  f is  the  mapping  indicated  by  (2,23).  This  mapping  is 
defined  in  the  subroutine  INDEX  which  is  used  by  BIGSET,  It  is  obvious 
that  f is  one-to-one,  and  therefore  f 1 exists,  f 1 is  defined  by 
taking  the  elements  of  the  solution  vector  v in  order  from  left  to 
right  and  placing  them  into  the  lower  triangular  matrix  L in  column  order. 

TTSTAR  is  then  called  to  find  the  lower  triangular  matrix  such 

that  M = T T*. 
r r r 

Now  that  T , exp(AtA),  and  T have  been  computed  and  the  vectors 
w^k^  for  k = 0,1, 2,  •••  can  be  constructed  from  the  w-sequence,  z(0)  can 
be  computed  by  (1.26)  and  the  sequence  of  vectors  z(kAt)  for  k = 1,2, •••  can 
be  computed  by  (1.29).  The  input  vector  SB  is  also  known;  therefore,  the 
specified  time  series,  x(t)  for  t = 0,  At,  2At,  •••,  can  be  computed  by 
(1.30).  This  process  is  carried  out  in  the  subroutine  SOLUTION.  This 
completes  the  construction  of  the  specified  time  series. 

The  card  input  to  GAUSSIAN  is  of  two  types,  problem  data  and 
flags.  The  sole  purpose  of  the  flags  is  to  control  the  flow  of  the  program. 
The  flag  NOSKIP  is  read  once  and  only  once  each  time  the  program  is 
compiled;  all  other  input  parameters  must  be  read  anew  for  each  time  series 
that  is  to  be  generated. 

The  following  is  the  list  of  problem  data: 

N n in  equation  (2.12).  N is  also  used  as  a flag 

in  that  GAUSSIAN  terminates  if  N = 0;  this  is 
to  be  used  for  normal  exit. 

M m in  equation  (2.12). 


sa( 1 ) , • ■ • , SA(N+l) 

V 

•••,  a^+1  in  equation  (2.12). 

sb( 1 ), • • ■ , sb(m+i) 

V 

•••,  b in  equation  (2,12). 

m+l 

DT 

At, 

time  interval  at  which  samples  of  the  time 

series  are  to  be  calculated. 
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IJOSKIP 


ISTOP 


I POWER  T 
I POWER  TR 


IPCWERB 

IPOWERC 


IFT 

IFTR 

IFB 


The  following  is  the  list  of  flags: 

the  number  of  end-of-files  (ie.,  previously 
computed  time  series)  which  must  be  skipped  on 
the  output  tape  before  the  first  time  series  is 
computed  and  then  buffered  out  onto  this  output 
tape . 

the  number  of  points  of  the  time  series  which 

are  to  be  buffered  out  onto  the  output  tape  in 

each  block  (except  for  the  last  block  which 

contains  the  remaining  points).  It  is  necessary 

that  ISTOP  S 5604. 

used  when  factoring  the  matrix  M. 

used  when  factoring  the  matrix  M^.  IPOWER  T 

and  IPOWER  TR  are  called  IPOWER  in  the  subroutine 

TTSTAR.  The  check,  "is  (m. . | < 1o"1POWER  ? /' 

is  made  before  m. . is  used  as  a divisor  in  TTSTAR 
11 

where  m. . are  diagonal  elements  of  M or  M . 
n r 

used  when  solving  the  set  of  equations  (l.2j). 

used  when  solving  the  set  of  equations  (1.26). 

IPOWERB  and  IPOWERC  are  called  IPOWER  in  the  sub- 

i -POWF'R  n 

routine  GAUSS  P.  The  check,  "is  | ^ | < 10  ?, 

is  made  before  a. . is  used  as  a divisor  in  GAUSSP 

ii 

where  a^  are  the  diagonal  elements  of  the  constant 

matrix  under  consideration. 

associated  with  IPOWERT. 

associated  with  IPOWER  TR. 

associated  with  IPOWERB. 
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I 

I 

1 

i 


IFC  associated  with  IPOWERC.  Let  IF  denote  any  of  the 

above  four  flags.  If  IF  = 1 and  a zero  divisor 
has  been  detected  by  its  IPOWER,  a diagnostic  will 
be  given,  but  the  calculation  of  this  time  series 
will  continue.  However.,  if  IF  = 0 and  a zero  divi- 
sor has  been  detected,  a diagnostic  will  be  given, 
but  the  program  will  go  to  the  next  time  series, 
print  the  number  of  points  which  are  in  the  fol- 
lowing block  of  time  series  points, 
do  not  print  the  number  of  points, 
print  the  computed  time  series  in  blocks  of 
ISTOP  points. 

do  not  print  the  time  series  points. 

The  data  deck  to  GAUSSIAN  must  be  stacked  as  follows: 

NO  SKIP  (Format  I 4) 


IF  N PRINT  1 

- 0, 

IF  X PRINT  = 1 

= 0, 


l 

l 

I 

j[ 

r 


» s 


i 


> 


N,  M ( Format  2l4) 

SA( 1 ) , • • • , SA(N+l)  (Format  4f20.6) 

SB(l),-"^  SB(M+l)  (Format  4F20.6) 

DT  (Format  F20.6) 

ISTOP,  IPO WERT,  IPOWERTR,  JPOWERB,  IPOWERC 
(Format  514) 

I FN PRINT,  IFXPRINT,  IFT,  I FT R,  IFB,  IFC 

(Format  6ll)  _y 

(Blank  card;  flags  normal  exit.) 

The  problem  data  are  printed  and  labelled.  The  intermediate 
matrices  and  vectors,  which  are  found  by  the  various  subroutines  are  printed 


Repeated 
for  each 

time  series 

to  be 

generated. 
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unlabelled  as  they  are  calculated. 

More  than  one  time  series  can  be  shaped  from  the  same  w-sequence 
each  time  GAUSSIAN  is  compiled.  This  input  tape  which  contains  the 
w-sequence  must  be  of  the  same  format  as  defined  for  the  output  tape  from 
CLLCNVRT.  Each  block  of  the  w-sequence  should  not  exceed  588  points. 
GAUSSIAN  terminates  when  and  only  when  either  N = 0 or  a parity  error  has 
been  detected  on  either  the  input  or  output  tape. 

The  output  tape  has  the  following  form: 

ISTOP 

X(l)  • • • • X( ISTOP) 


ISTOP 

X(l)  • • • • X( ISTOP) 
IPRINT 

X(l)  * * * * X( IPRINT) 
end-of-file 

ISTOP 

X(l)  • • • • X( ISTOP) 


IPRINT 

X(l)  • • • • X( IPRINT) 


end-of-file. 


III.  "WHITE  SEQUENCE, " w-SEQUENCE,  AND  EXAMPLE  TIME  SERIES 

A spectral  analysis  program  similar  to  the  program  CPECT,  as 

15 

found  in  a report  by  Ellis  and  Boston,  is  used  to  find  the  autocorre- 
lation function  and  power  spectral  density  function  of  the  "white  sequence" 
generated  hy  ISLASHO  and  four  example  time  series  generated  by  GAUSSIAN. 

The  autocorrelation  function  is  defined  by  (1.3)  and  the  power  spectral 
density  function  is  defined  by  (1.4).  In  SPECT  the  time  t in  (1.3)  and 
(1.4)  is  taken  at  regular  intervals  At,  and  therefore 

r = IAt,  1=0, 1,2,  •••,  (3.1) 

where  1 is  called  the  lag  number.  For  practical  applications  on  the  com- 
puter, t can  only  take  on  finite  values  and  the  sequences  considered  will 
only  have  a finite  number  of  points;  therefore  a finite  number  of  lags 
will  be  taken  in  SPECT.  The  maximum  lag  number  will  generally  be  about 
ten  percent  of  the  number  of  points  in  the  given  sequence. 

A.  "White  Sequence" 

The  condition  (l.l4)  states  that  the  autocorrelation  function 

of  a "white  sequence"  must  be  zero  for  all  but  the  zeroth  lag.  The 

"white  sequence"  under  consideration  is  le' ),  1 £ n § 4266.  It  is  also 

stated  in  Chapter  I,  Section  B,  that  the  power  spectral  density  function 

of  v must  be  constant.  The  plots  of  the  autocorrelation  function  a.nd 
n r 

the  power  spectral  density  function  of  vn  as  computed  by  SPECT  for  tim® 
interval  At  = . 1 sec.  and  1 = 200  lags  are  shown  in  Figures  1 and  2, 
respectively.  It  can  be  seen  from  these  two  plots  that  the  sequence 
(enJ  is  a good  numerical  approximation  to  a "white  sequence." 


30 


B.  w-Sequence 

A short  program,  which  is  not  listed,  was  written  to  find  the 
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mean  and  the  variance  of  the  w-sequence.  Let  the  mean  be  m = E(w^)  and 
the  variance  be  defined  as 


* i L (v 


(3.2) 


i-1 

where  w^,  1 S i S H,  are  the  N points  of  the  w-sequence.  It  Is  stated  in 

Chapter  I that  the  w-sequence,  as  constructed  by  (1.17),  should  have  mean 

0 and  variance  1.  The  numerical  w-sequence  constructed  from  { e n j by 

2 

CLLCNVRT  was  found  to  have  m * - .029  and  o * .9577. 

C . Four  Example  Time  Series 

Four  example  problems  are  now  to  be  considered.  The  following 
three  steps  are  taken  in  the  consideration  of  these  examples: 

(1)  Given  the  constant  coefficients  of  the  polynomials 

P(  iui)  and  Q(iu>)  of  the  theoretical  power  spectral  density 
function  V(u>)  as  given  by  (2.12)  and  a time  interval 
At,  GAUSSIAN  is  used  to  shape  a time  series  with  this 
power  spectral  density  function. 

(2)  A log-log  plot  of  the  theoretical  power  spectral  density 
function  as  a function  of  frequency  f,  is  made  where 

(i)  . 

? — and  0)  is  the  angular  velocity  given  in  radians/sec. 

(5)  The  actual  autocorrelation  function  and  power  spectral 
density  function  are  then  calculated  by  SPECT;  the  auto- 
correlation function  is  plotted  linearly  against  time 
while  the  power  spectral  density  function  is  plotted 
against  frequency  on  a log-log  plot. 

Only  the  power  spectral  density  function  is  considered  in  the  first 
three  examples;  however,  both  the  autocorrelation  and  the  power  spectral 


I 

i 
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density  functions  for  example  four  are  given  by  Richie^  and  will  therefore 
be  considered. 

In  all  cases  the  theoretical  and  actual  log-log  plots  of  the 

power  spectral  density  functions  were  found  to  have  approximately  the  same 

shape;  however,  for  each  example,  the  two  curves  were  found  to  differ  by 

some  arbitrary  constant.  If  the  magnitude  of  the  points  of  the  power 

spectral  density  function  is  specified  and  not  just  the  shape, 

the  specified  time  series  can  be  obtained  by  multiplying  each  point  of  the 

original  time  series  by  the  square  root  of  the  constant  whose  logarithm 

is  the  difference  between  the  theoretical  and  the  actual  curves. 

It  can  be  seen  by  (1.3)  that  if  each  point  of  the  time  series 

x(t)  is  multipled  by  a scale  factor  s,  then  for  every  t,  the  autocorre- 

2 2 

lation  function  is  multiplied  by  s . Since  the  constant  s can  be  removed 

from  the  integral  in  (1.4),  the  power  spectral  density  function  is  also 
2 

multiplied  by  s . 

Let  t^  and  a^  be  corresponding  points  on  the  theoretical  and 
actual  power  spectral  density  curves,  respectively,  and  if  there  exists 

an  s > 0 such  that  for  every  i where  t^  and  a^  are  defined, 

2 2 

log  s = log  t^  - log  a^,  then  t^  = s a^.  The  converse  is  also  true. 
Therefore  the  specified  magnitude  of  the  power  spectral  density  function 
can  be  obtained  by  multiplying  each  point  of  the  time  series  by  the  scale 
factor  s.  This  scale  factor  may  have  physical  units. 

In  Figures  3 through  6 the  crooked  continuous  curve  is  the 
actual  power  spectral  density  function  as  computed  by  SPECT  and  the 
smooth  dotted  curve  is  a multiple  of  the  theoretical  curve  which  gives 
an  approximate  fit  to  the  points  of  the  actual  power  spectral  density 
function. 
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The  theoretical  power  spectral  density  function  of  example  one 
is  given  to  be 

V (oj)  = , (3-3) 

<u  -60/  + 25 

which  has  zeros  at  the  points  (0,l/3)  and  (0,-l/3)  and  poles  at  the  points 
r,  r,  - r,  and  -r  where  r = (2,1).  Therefore  for  all  real  a>,  the  condition 
0 < V^(u>)  < 00  holds.  Also  for  real  ou,  the  conditions  V^(ur)  = V^( -oj)  and 
V^o))-*  ± °°  hold;  therefore  the  three  conditions,  as  stated  in  (1.8),  exist 
which  insure  that  V^(u>)  can  be  represented  in  the  form  of  (2.12).  Specif- 
ically, V^cn)  can  be  represented  as 


vL(«)  = 


3(ia>)  + 1 


( icu)  + 2(  io>)  + 5 


(3-4) 


where  1+25+5  has  zeros  at  (-1,2)  and  (-1,-2).  Both  of  these  zeros 
lie  in  the  halfplane  Re  5 < 0. 

The  theoretical  power  spectral  density  function  of  example  two 
is  given  to  be 

k 2 s 

...  i ~z  Cr\ 

(3.5) 


. , / \ cd4  - 70o>2  + 1369 
v2(  “)  = ~Z — 774 


oj  + 14<jo  + 4900  + 36 

The  numerator  of  V,,(a>)  has  zeros  at  the  points  r,  r,  -r,  and  -r  where 

r = (6,l),  and  the  denominator  of  Vg(ao)  is  greater  than  zero  for  all  real 

00;  therefore,  for  all  real  00,  the  condition  0 < Vg(cjo)  < « holds.  Also  the 

other  two  conditions  of  (1.8)  hold;  therefore  can  be  represented,  as 


indicated  in  (2.12),  by 


V2(m)  = 


(ia>)  + 2(ia))  + 37 


(ioo)^  + 6(iui)2  + ll(io))  + 6 
2 2 ' ' 
where  5^  + 65  + 115  + 6 has  zeros  at  (-1,0),  (-2,0),  and  (-3,0).  All 


(3.6) 


these  zeros  lie  in  the  halfplane  Re  5 < 0. 


The  theoretical  power  spectral  density  function  of  example  three 
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is  given  to  be 


„ , \ 40Qas  + 1 

V (a>)  = -7 u 2 

as  + l4as  + 4903  + 36 


(5-7) 


V^(a>)  has  zeros  at  the  points  (0,  and  (0,  - ^).  The  denominator  of 
this  example  is  the  same  as  in  example  two.  It  can  be  seen,  as  in  examples 


one  and  two,  that  V^(us)  can  be  represented  by 


Vx(«d) 


20  ( ias)  + 1 


(ito)^  + 6(ios)^  + ll(ia>)  + 6 


0.8) 


Example  four  is  an  example  problem  which  appears  as  channel  4 in 
17 

a report  by  Richie.  The  theoretical  power  spectral  density  function  is 


given  to  be 


2 / , 2 2, 

v f(o)  - 4Ak-n M + (k  + c ) 

V4v“;  - 2 2x  2 ~2  2x2 

o>  + 2(k  -c  )o3  + (k  + c ) 


(5.9) 


where  A = 2149.2,  k = .658,  and  c = 2.71.  Since  A,  c,  and  k are  positive 
and  k > c,  the  conditions  (1.8)  hold  for  real  us.  Therefore,  V,(os)  can  be 


represented  as 


V.  (as)  = 


b1(ias)  + b2 


a1(ias)  + a2(ios)  + a^ 


(3-10) 


where  ^ * 75.210999497,  t>2  ~ 53.521755647,8  = 1.0,  ag  = 1.316,  and. 
a,  = .506405. 


Both  zeros  of  a^l  + a.^  + a^  lie  in  the  halfplane  Re£  < 0 since 
a,  > 0,  a >0,  and  a > 0.  This  fact  can  be  shown  by  letting  x = Re£  and 

C.  J 

y = Im£  and  then  assuming  that  there  exists  a 5 with  x £ 0 such  that 
2 

a^S  + = 0.  This  complex  equation  is  equivalent  to  the  following 

two  real  equations: 


iL(x  -y  ) + a2x+a^ 


(3.11) 


2a ^xy  + a2y  = 0 


(3.12) 


MM 


HM5  e . ITT  ; 

••  .If  ~ 
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p 2 

By  (3.12),  y = 0 since  x 2 0.  Then  by  (3.1l),  a^(x  -y  ) + dgX+a^O 

2 

which  is  a contradiction;  therefore  both  zeros  of  a,  i + + ^-n  the 

halfplane  Re|  < 0. 

The  actual  power  spectral  density  function  of  the  time  series 
generated  by  GAUSSIAN  for  this  example  was  found  to  differ  from  the  theoret- 
ical power  spectral  density  function  only  by  a multiplicative  constant; 

. 18 

however,  the  theoretical  autocorrelation  function  given  by  Richie  and  the 
actual  autocorrelation  function  differed  greatly.  This  illustrates  the  fact, 
as  pointed  out  previously,  that  small  perturbations  in  V(co)  may  result  in 
great  changes  of  R(t).  These  two  autocorrelation  functions  are  shown  in 
Figure  7,  and  the  power  spectral  density  functions  are  shown  in  Figure  6. 

The  degree  and  coefficients  of  the  polynomials  P and  Q,  which  are 
used  as  input  parameters  to  GAUSSIAN,  can  be  taker,  from  equations  (3.M, 

(3.6 ),  (3.8),  and  (3.IO).  These  parameters,  along  with  the  time  intervals 
used  in  generating  a time  series  for  each  of  these  examples,  are  also  given 
in  Table  1. 
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77  PRINT  7b 

7«  formai  (?<sh  Parity  frrofi  on  mag  Tapf> 

<79  CONTINUE 
F NO 

3f>0()  Fortran  otagnosiic  rfsolts  - for 


ISl.  A SRI) 


NO  errors 

COMPASS. L.X 


COMHASS- V 

< 2>. 

1 ) 

WANOOM 

oonoo 

01077777 

0| 

n 

77777 

0 

RANDOM 

ENTRY 

UJP 

RANDOM 

»• 

UUQ01 

14600000 

14 

1 

ooooo 

2 

ENA 

0 

0000? 

40000007 

4 0 

0 

onono7 

0 

ST  A 

CARRY 

0000  3 

14177776 

14 

0 

777  76 

1 

ENT 

-1  . 1 

00004 

20000005 

20 

0 

0000(16 

0 

LOA 

L 

00006 

164  77776 

16 

1 

777  76 

0 

INAiS 

-1 

iionoh 

44000007 

44 

0 

POOOO 7 

0 

Ska 

•♦1 

00007 

10177777 

10 

fl 

77777 

1 

INUt X 1 

ISt 

00010 

01000012 

01 

0 

000012 

0 

U Jo 

• ♦? 

0001  1 

01000022 

01 

0 

00002? 

0 

UJO 

F1NISH1 

0001? 

20100000 

20 

0 

cnonno 

1 

LOA 

1 1 

U0013 

50000002 

60 

0 

000002 

0 

Ml  IA 

17  1 

00014 

30000007 

30 

0 

OOOOO  7 

0 

ADA 

CARRY 

oqq35 

13077747 

33 

0 

77747 

a 

Shad 

-24 

00016 

61000006 

61 

0 

000006 

0 

OVA 

THOUSAND 

00017 

41127340 

41 

0 

C27340 

1 

STO 

IY,  1 

000?0 

40000007 

40 

0 

noono  7 

0 

STA 

carry 

00021 

01000007 

01 

0 

000007 

0 

UJP 

INUEX1 

000?2 

2000000S 

20 

0 

000005 

0 

F TNT5H1 

LOA 

L 

UP  023 

53500000 

53 

1 

noouu 

1 

IAI 

1 

00024 

20000007 

20 

0 

000007 

0 

LOA 

CARRY 

0002b 

40127340 

40 

0 

C27340 

l 

STA 

IY.l 

00026 

1 4600000 

14 

1 

ooooo 

2 

ENA 

0 

00027 

40000007 

40 

0 

D00007 

0 

STA 

CARRY 

000.30 

141 77776 

14 

0 

77776 

1 

ENT 

-1.1 

niwni 

200UQ005 

20 

0 

onoo ob 

0 

LDA 

L 

00032 

15477776 

15 

1 

77776 

0 

INA  .S 

-1 

00Q33 

44000034 

44 

0 

P00034 

0 

Ska 

*♦1 

00034 

10177777 

in 

0 

77777 

1 

INDEX? 

1ST 

*••1 

00035 

01000037 

01 

0 

*00037 

0 

UJP 

•»? 

00036 

01000060 

01 

0 

000050 

0 

UJP 

F 1 NlSM? 

OOQ3J 

203  QA1)  OP 

20 

0 

CQQOOO 

3 

LOA 

IX. 3 

00040 

60000003 

so 

0 

000003 

0 

MU  A 

IX? 

00041 

30000007 

30 

0 

1)00007 

0 

AOA 

CARRY 

00042 

30127341 

30 

0 

027341 

1 

AOA 

l V. 1 , 1 

00043 

13077747 

13 

0 

77747 

0 

SHAU 

-24 

00044 

61000006 

61 

0 

000006 

0 

OVA 

thousand 

(10045 

43327341 

41 

0 

C27341 

J 

sia 

IY.1.3 

00046 

40000007 

40 

0 

000007 

0 

STA 

carry 

00047 

01000034 

ol 

0 

O00034 

0 

UJP 

INDEX 2 

OOOSO 

20000006 

20 

0 

1)00006 

0 

F 1 N I SH? 

LOA 

l. 

00051 

16600001 

16 

1 

00001 

2 

INA 

t 

00052 

53500000 

63 

1 

nnooo 

1 

TAT 

1 

axui53 

20000007 

20 

0 

000007 

u 

LOA 

CARRY 

00054 

40127340 

40 

0 

C?  7 340 

1 

STA 

IY.l 

00055 

14600000 

1 4 

1 

ooooo 

2 

ENA 

0 

00056 

40000007 

40 

0 

000007 

0 

STA 

CARRY 

00057 

141/7776 

14 

0 

777  76 

1 

ENI 

-1.1 

00060 

20000006 

20 

0 

000006 

0 

LOA 

L 

00061 

15477776 

15 

l 

nut> 

0 

LNA.S 

-3 

00062 

44000063 

44 

0 

000063 

0 

Ska 

•♦1 

00063 

101 77777 

10 

0 

77777 

1 

IN0FX3 

IST 

*•.1 

00064 

01000066 

01 

0 

000066 

0 

UJP 

*♦? 

02/0 


( 


47 


COMF'ASS- 1? 

(?.  1 ) 

HIM 

OOObS 

0 1 Oil 0 0 7 7 

01 

0 

O00077 

0 

UJP 

FINISH)  )/bb 

oonoo 

?<>1  HOOOO 

2(1 

0 

r n o o o o 

1 

L lift 

1 > . 1 

U(IOb7 

5:1000004 

0 

000004 

0 

Ml  1 A 

1 * 3 

oon  7o 

3000000 7 

io 

0 

000007 

0 

ADA 

canny 

onn  7 1 

3ii  1 ? 7 34? 

0 

C 2734? 

l 

AD* 

lv*2.1 

000  7? 

1 <0  7 7 74  7 

1 1 

n 

7774  7 

0 

Shad 

-24 

00073 

S 1 onoonh 

M 

0 

00000b 

0 

DV* 

thousand 

000  74 

41 1 27  34? 

41 

0 

C 3734? 

1 

STn 

IV*?, 1 

00075 

4000000  7 

4 0 

n 

1)0000  7 

0 

ST* 

CAPPY 

000  7b 

0 1 000i)r>3 

0 1 

0 

00 00b 3 

0 

UJP 

INDEX  3 

00077 

2000000 7 

/Ml 

0 

1)0000  7 

0 

F 1NISH3 

LIIA 

canny 

ooi  on 

0 <1  00  1 05 

0 

Pool  05 

1 

A Z 1 , NF 

7 F wo 

0010) 

20000005 

20 

0 

1)00005 

0 

LDA 

1. 

0010? 

lbbllOOO? 

is 

1 

00002 

2 

INA 

? 

00103 

40000005 

40 

0 

nnooos 

0 

S1A 

L 

00104 

0 1 000 1 1 4 

0 1 

0 

PoOl  1 4 

0 

UJP 

change 

UOIOS 

20000005 

20 

0 

000005 

0 

/FPO 

LI)A 

1. 

0010b 

lbbOOOO? 

IS 

1 

00002 

2 

IMA 

2 

00107 

S 3500000 

S3 

1 

00000 

1 

3 At 

1 

OOI  10 

15000001 

is 

1 

00001 

? 

INA 

1 

00111 

40000005 

4 0 

0 

000005 

0 

STA 

1. 

0011? 

20000007 

?n 

0 

000007 

0 

LDA 

caphy 

00113 

40  12  7 340 

4 0 

0 

C27340 

l 

STA 

TV.l 

00114 

|4l 7777b 

14 

0 

777  7b 

l 

CHANGF 

F N T 

-1 .1 

001  IS 

20000005 

?0 

0 

1)00005 

0 

LDA 

L 

OOI  lb 

1547777b 

is 

1 

7777b 

0 

1 N A , S 

-1 

00117 

44000120 

44 

G 

pool  20 

0 

SWA 

1 

001?0 

101 77777 

10 

0 

777  7 7 

1 

INDF  X 4 

IS? 

»*.l 

00121 

01000123 

01 

0 

POO  1 23 

0 

UJP 

»♦? 

001?? 

0 1 0001 ?h 

0 1 

0 

POO  1 2b 

0 

II JP 

F IMISH4 

00123 

20127340 

?0 

0 

C27340 

1 

Ln* 

1 Y * 1 

00  1 ?4 

40100000 

4 0 

0 

rooooo 

1 

STA 

IX, 1 

00125 

01000121) 

01 

0 

POU120 

0 

UJP 

1NUEX4 

001?b 

2000000 l 

20 

0 

OOOOO 1 

0 

F INISH4 

LDA 

1DFC 

00127 

30000000 

10 

n 

OOOOOO 

0 

ADA 

NODEC 

00  1 30 

40000001 

40 

n 

OOOOOl 

0 

ST  A 

idfc 

00131 

0 1 000000 

01 

0 

PilOOOl) 

u 

UJP 

DATA 

panoom 

00000 

NOUEC 

HSS 

1 

00001 

If)  FT 

HSS 

1 

0000? 

1*1 

HSS 

l 

0000  3 

IK? 

HSS 

1 

00004 

1X3 

HSS 

1 

oooos 

L 

HSS 

1 

00005 

I HO') 5 AND 

ass 

1 

0 0 0 0 7 

CAWOy 

HSS 

COMMON 

1 

00000 

lx 

HSS 

1 2000 

27340 

ir 

HSS 

12000 

END 

NUMHFP  OF  l INFS  with  DIAGNOSTICS 


1*8 


*?00  FOKTWMv!  (?.  I ) 


PP(n  P,.  > (II  f I H/W  I 
I I-  HP  •»(»  I (II) 

? kip.-.ai  i^km  paoiiv  ikmim  >>„  (nitPul  Tap*- 

( f-'lP-AI  (^/M  PA4I  If  IWPlIH  n.  I'MPllt  iapci 
4 F HP  "A  I (/////I 
PI-  All  la  II-  pm  | |VI  I 

CALI.  UINVI-MI  III  UP  I N T a I PAMlIr  la  I 
IF  < I PAM  IT  V I .11.  I)  I .)  a II 
1(1  PM  In  | ( 

PMTi-I  < 

(ad  10  -)P 

11  IF  (1  PAMllr  (I  ,F0.  1)  l?a  MM 
1?  PM  IN  I ,» 

PMTM  4 
MM  CONI  1 ytll 
t No 

1^0(1  piIpTMAj  'I  I AfaNOS  | IC  pFS'ILTS  - FOP 


Ol/ni/FiF, 


I 


PArillY  O) 


CLI-CNVM  T 


NO  FRRORS 


lt9 


[ 

I 

I. 

1 

I 


1 

1 

l 

[' 

( 


mwihhi  (?.  i ) n vovhh 


SOHwuoi  | COHVFWT(  If  HIM  I t.  | haw  1 T Y I,  T WAWITY  <)) 
OlMhNSlU'l  X(bO')> 

WfAl  LOOP 

1 POP  'A  I (I  | it) 

2 FOWmA  I (tit/  > . 1 0 t 

3 POH'IAI  (//) 

S FftH'iAl  ( I H|t 

KtWlNJ  I 
PP  w 1 Pit i s 
I PAH  i I Y 1 = 0 
I PAHIIY  (I  = II 
PWlNT  S 

20  HOpPtW  1M  (1,1)  ( 1 HOF.  IHOP) 

21  GO  10  (21.  22,  23.  77),  UFlIISTFlb) 

23  t NO  f ILF  1 
PWTimT  5 
WFTUPn 

77  I PAHITY  I = I 
wFTown 

22  l BuP  = 1 bUF 

HUFFtw  1M  (S. 1 ) ( * ( l ) ,X  ( I Hof)) 

24  GO  10  (24.  2b.  23,  77),  UNlTSTF(b) 

25  I HOF  2 = I HOF/ 3 
PI  = J.1415W2653S 
TPT  = 2 • II*  PI 

00  2B  1=1,  I tUJF  2 
12=1*2 
I2Ml  = 12-1 

X 1 = Sow  1 P ( -2. 0*|.  DGF  ( X ( 1 2M  | ) ) ) 

X2=  Tw 1*  X ( 12) 

X ( t 2m 1 ) = x 1 «roSF ( x A) 

28  XU2)=X1“SINF  (X2) 

HOFFFw  ()|| I (1,1)  ( 1 HOP.  1 HOF) 

31  GO  10  (31.  32,  32,  7H>.  UN(TSTF(11 

78  l PAwnr  o = i 
pf  town 

32  I HOF  = 1 HOF 

HOFFFw  our  0,1)  (X(l),  X(I  PDF)) 

33  GO  To  (33,  (4.  14.  78),  IMtTSTFd) 

34  I F ( [F  PW I NT  .Ft}.  1)  4il,  20 
40  PW I NT  I,  1 HOF 

PWInT  2.  U(I>,  1=1,  I HOF) 

PW  1 FIT  3 
GO  TO  20 
FNI I 

3200  FOWTWAN  OlAGNOSflC  WFSoLTS  - FOP  CONVPwT 


NO  EWWOPS 
F0UIP.1=MT 
LOAD, 56 
PUN 


I 


l 


3200  Fortran  (2.11  3 1 /'i  3/66 

PROGRAM  GAUSSIAN 

COMMON  SA  ( 1 3)  « SR ( 1 2 ) • C(7H.  7R)  . «(7«).  T(l?.  12).  F (1  ?.  )?), 

I A 1 1 2.  12).  6112.  12) 

1 FORmAT  ISUI 

2 FORMAT  !4F2n.G) 

3 FORMAT  (///,  aSH  TMF  M MATRIX  IS  SINGULAR) 

4 FORMAT  (///) 

5 FORMA!  (SF?'.10) 

S FORMAT  C///.  2SH  RIGStT  MATRIX  IS  SINGULAR) 

7 FORMAT  (///.  m7H  PAHllY  FRROR  ON  INPUT  TAPF) 

R FORMAT  (///.  2HM  PARITY  ERROR  ON  OUTPUT  TAPt) 

R FORmATT///,  2JH  T TSTAK  OOFS  NOT  fXTSYT 

10  FORMAT!///,  2SH  TR  TRSIAR  DOES  NO  1 EXIST) 

11  FORMAT  (hill 

12  FORMAI  (14H  DETERMINANT  = .E20 . I U . // ) 

100  FORMAT  I 1 H 1 ) 

SS  FORMAT ()oh  CONSTANTS  OF  DENOMINATOR  POLYNOMIAL  ) 

47  format  ( I4H  CONSTANTS  OF  NUMFRATor  roi  YnOmIaL  ) 

M3  FORMAT  («Fl  N m ) 

SR  FORMAT  (SX.  \3HTTmF  INTERVAL  ) 

70  FORMAI  (10X.  SK(SFC)  ) 

RF  w t no  > 

RtAU  1.  NO  SMP 
HO  ..  | 1 = 1 , NO  SNIP 

S3  PUFF  ER  IN  < ? , I ) ( I SK I P . ISMP) 
s?  GO  TO(S2.  S3.  S'.  S4).  UNI1STFT2) 

6m  PR  I of  6S 

SS  F0RmAT(4"H  PARITY  FPROK  On  mai,  | apF  WMlLF  SkTPpINi,  ) 

1.0  TO  Ft 
S!  ('ONflNUF 
SO  RtAU  1.  N*  M 

IF  IN  .Ml,  ) RR.  71 

71  NP 1 = N«  1 

MP 1 = M*1 

PHI  iT  100 
PRT'.T  -m 
PRINT  1 . N.  M . 

REAO  ?,  (SA ( I ) . T = l . NP  1 ) 

REAU  2.  (Sail).  1=1.  MPT) 

READ  Of 

REAO  I.  ISTOP.  IPoWtHI.  IPOWERTR,  IPO-tRB.  IPowFHC 
REAO  i).  IF  N PR  T Ml  « IF  X PRINT.  (FT.  TFTP.  I F i- . TFC 
PRTuT  A 
PR  1 it  T SS 

PRINT  5.  ( SA ( I ) . 1=1.  NP1) 

RR 1 nT  « 

PRInT  t-7 

PRInT  S,  (SH ( l I . 1=1 . MP1) 

PRINT  4 
PPTnT  SR 
PRIafI  70 
PRINT  S.  DT 
PRInT  4 
CALL  FIND  H(N) 

DO  SO  1=1.  N 

SO  PRINT  S.  (C(I.J).  J*l«  NP1) 

PH  In  I *. 

CALL  GAUSS  P (N.  KS.  IPOwFHH,  DET) 

PRINT  12.  DET 


50 


f 


51 


It-  IKK)  21  1 20 
?\  PRINT  3 

If  < If  H .fU.  I ) ?0.  60 
?0  CALL  FIND  CM(N) 

PRINT  St  Ulllt  I a 1 « N) 

PRInT  a 
DO  SI  1=1 t N 

SI  PH  TNT  St  (C ( I t J>  t J= I t N) 

CALL  I T STAH  ( Mi  KKt  IPOwFRT) 

IF  (KM  Alt  A'l 
A)  PRIM  <i 

IF  < IF  I ,FO.  1|  A 0 1 60 
An  CALL  MAKf  AIN) 

PH  I NT  A 
DO  52  1= l t M 

SP  PH  TNT  St  ( I ( 1 t J) t J3 I , N) 

PRINT  a 
DO  S3  I = | t N 

S3  PRINT  St  ( A I I t J ) t J=lt  N) 

CALL  tXP  A 13  T ( N , DTI 
PRINT  a 
DO  SA  I = \ t N 

5a  PRINT  5.  IEU.J).  j:|t  N) 

CALL  H I GSfc  TIN) 

NN  3 lN»(N*l)>/P 

NNP13  NNt 1 
PRINT  a 
DO  S5  I = t , NN 

55  PHIlJ-T  St  ICIItJlt  J=lt  NNP1) 

CALL  GAUSS  P I NN  t KKt  IPOwFRCt  OF  I ) 

PRINT  a 
PRINT  lif.  OFT 
IF  IKK)  22t  P.3 
??  PH  I n T 6 

IF  HFC  .tli.  1)  ?3t  fai) 

23  I 3(1 

PRINT  a 

PRINT  St  (XII) t 1=1 t NN 
DO  2 A 1= I t N 
DO  a A J=lt  N 
L*L*i 

( ( I t J) = X IL) 

?A  C(JtI)=  X (L ) 

PRINT  A 
00  56  I*  1 9 N 

56  PRINT  St  (C I I t J)  t J= I t N) 

DO  25  1=1.  N 

DO  25  J=1 t N 
PS  A 1 1 , J) = 1 <I,J) 

CALL  I I SlnRlNt  KKt  IPOwfRTR) 

IF  IKK)  A3 1 A2 
A3  PRINT  in 

IF  I IF  TP  .E'J.  ll  A?t  6 U 
A?  HtWlNl)  I 
PRINT  A 
DO  ST  I * I t N 

57  PHINT  S.  IT  I It  J) . J=l . N) 
print  inn 

CALL  SOLU T ION (NtRtlPAHlTYItlPARl T TO. IS TOP t IFNPR1NT .IFXPHINI ) 
IFllPARITVI)  26 1 <*8 
26  PRINT  ? 

GO  TO  HR 


I 


I 

I 

9R  IF  ( IPARI  TY1»)  ? T.  f.0 
PRTnT  H 
99  PR  I >9 T 100 
F.NI 

3^00  FORTRAN  1 AGNOSTIC  RESULTS  - Fr\R 


UNDEFINED  SIMPLE  VAR1AHLFS 
IblUF 


♦ 


GAUSS  I AN 


i?on  hwiuiN 


(Z.  1 ) 


SUhmOU 1 1 Nt  F INI)  ' iM) 

COMMON  SA 11 1) «Sm 1 1?) .M I 7H.79) 

INTFGFP  0,  <JM| 

NM)=  N-l  * NP  | = N ♦ ) 

III)  I 1 = 1.  NM 1 
1 F>(I.  NP  11=  0.0 
MIN.  NP1)=  •>.*, 

00  z I=|,  N 
1)0  Z J=l,  N 
Z HIT.  J)=  0.0 

no  ^ k= i . n 

K M 1 = K—  1 % KL  = KMl  /<>.  II  ♦ . S 1 *1 

MJ  = (N«NMl)/2*l  * A3=  <-1.0)««Hrtt 

III)  I ll  = KL.  Ml 
UM)  = I)-)  i N4=  N— 2*UM  | *K 
3 H (K .U) = A3* |-1 .0) «»OM)«bA (M.) 
tN.l 

.IZOll  FOHTHAN  liTAOMlbT  IC  PtbilLlb  - FoP 


F |Ni)F> 


NO  FPMOPS 


r' 


3?oo  t-oJTW'N  ap.ii 

S'lW.OilIlNt  flNllC*l(N) 

COMMON  SA  < 111  .son  ?»  «CM<  7«,  7P>  . SM  < 7H) 
l>0  l 1=1*  N 
1)0  I J«l«  N 
l CM(t*J|*  u.lt 
I >0  ? III.  N 
Ar  —l.i) 

DO  > Jr  I.  N.  >< 

Ar  -A  4 Kr  (J.J)/? 

CM(I.J)*  A * SM ( K ) 

? CM ( J» l ) * CM ( I . J) 

FNm 

3S00  F OPTPAN  DIAGNOSTIC  MFSULfS  - F OR  F I NOT  M 

NO  ERPOPS 

I 

I 

1 
I 

[ 

I 

! . i 

i 

; it 

: !l 


L_  L_ji- 


54 

31  /f'l/hfs 


55 


.1 


l 


I 

1 


j 


i 

I 


I 


1 

1 

l 

1 

I 

r 

i 

i 

i 

\ 

i 

i 


3?on  fortran  <?.ii 


1 1 /<>3/hF 


SURHOUTINE  T r STAR  (N,  KK  » IPOsFRT) 

COMMON  SA  ( ) 1)  ,SR(  I?)  «CM(7H«7'»)  ,SM(  78)  .T  (1?.  I?) 


C THIS  SURHOUTINE  lb  DESIGNED  TO  FACTOR  i>  POSITIVE  OF FINITE, 

C REAL,  SYMMETRIC  MATRIX,  Cm,  InTO  <.  LO -Fh-T w I ANGWL AR  MATRIX,  T. 

C V 1 Tm  POSITIVE  El FmFNTs  ON  THE  MAIN  DIAGONAL,  SUCH  THAT  CMsf  T». 


C KKsO  IF  THE  TRIANGULAR  MATRIX  SOLoTINN  EXISTS,  IF.  «LL  THE 

C ELEMENTS  ON  The  main  DIAGONAL  of  THE  rM  MaTR|»  ARE  POSITIVE. 

C KKa  I IF  THF  SOLUTION  IIOES  NO  I FXtsT. 


KKs  . 

1 £N  POtfFfc  = 1 0 • n « « (-IPUWFRT) 

00  H 1=1,  N 

IF  ICM(l.I)  ,LT.  TEN  POUEM)  t . s 
T KKsl 
GO  TO  * 

F CONTINUE 
A NM 1 = N-l 

00  S 1 = 1 , NM  1 
IPT*  I • I 
DO  H J=IP1 , N 
S T ( I . J) = 0.0 

T (1 ,1 )=  SOHTF (CM ( l , 1 1 ) 

UO  1 1=2,  N 

1 T ( T , 1 1 = CM ( I , | ) / T 
UO  ? Jx?«  N 
SUM=  0.0  f JM1=  J-l 
DO  3 K = 1 . JM 1 
T SUMS  SUM,  T ( J.K ) *T ( J,K ) 

T (J.J)sSuHTF  U.M1  J.Jl-SUM) 

JP1=  J*1 

DO  2 IsJRJ.  N 

SU«=0.n 

1)0  A Ksl.  JM  1 
A SUM=  SUM, T(T«K)*T(J,K) 

? T (I , J) * (CM( I , J) -SUM) /T ( J, J) 

END 


3200  FORTRAN  DIAGNOSTIC  RESULTS  - For  ttstar 


NO  ERRORS 
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3200 


(?.  ) ) 


3 1/03/66 


SUHMOOT1NF  MANF  AIN) 

COMMON  SA  ( 1 1)  ,SM  ( 1 ?)  ,C  ( 7H,  79)  . * < 7«>  , T I 1 ?.  )*>)  ,F  I I ?,  1?)  , A ( 1 1?) 
NM1=  '1- 1 
NP?=  n*,> 

DO  | J = l,  MM) 

DO  I J- 1 . N 
1 A(1.J>=  0.0 
DO  i 1=1.  NM1 
1P)=  1*1 
1 A (I  ,IP1 ) = 1 .<• 

DO  A 1=1.  N 

Nl*  NP^-| 

A A<N,I)=  -SA(Nl) 

EN'  i 


3/>00  FOMT9AN  DIAGNOSTIC  WF  SDL  T S - Fn»  maKFA 


NO  F RHOBS 
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3200  FORTRAN  <?.l) 


19/OA/hh 


SUHWOUUNt  FXP  A OTIN.DT) 

COMMON  SA I 1 1) .5HII?I *AN( 12.12) ♦ «! 12.1?) .DUMMY (SH7 A ) . X ( 78 ) , T ( | ? , 1 ? ) 
l «i(12<l?)<  A I 1 2 , 1 2 > 

NM 1 = N-l 
FACN=  1.0 
DO  1 I a 1 . N 
DO  1 J=l,  N 
AN  ( 1 . J)  = A t 1 . J) 

1 fc ( I . J) = 0.0 
DO  2 1=1.  N 

2 fc  1 1 » I ) * 1.0 
EPS CMK  = 0.0 
DO  7 1=1.  N 
CHfcCb=  AHSF I AIN. 1 1 ) 
lFICHtCb  .faf.  FPSCHK)  a, 7 

8 fc PSCHis  = CHfcCK 
7 CONllNUt 

fcPSCHb  = fcPSCHK«l .Ofc-12 
S TF  N=  OT/FACN 
00  12  1=1.  N 
DO  12  J= 1 , N 
12  HI  I . J) = AN! ! , J) « TFN 
DO  2H  J= 1 , N 
DO  2H  1=1.  N 

IF  I AHSF (HI  1 . J) ) ,LT.  tPSCHM  PH,  9 
2H  CONI lNOfc 
(.0  TO  99 

9 DO  b 1=1.  N 
DO  b J= 1 . N 

b F.  1 1 . J)  = F 1 1.  J)  «H(  I • J) 

DO  11  1=1,  NM 1 

1P1  = 1*1 
DO  II  J= | , N 
II  AN(l.J)  = HC1P1.JI 
HNN  = H(N.N) 

AN  I N , 1 ) = HNN*A I N, 1 ) 

DO  10  1=1,  NM1 
1P1  = I ♦ 1 

10  AN(N.lPl)  = HIN. 1 1 ♦HNN*AIN, 1PI > 

FACN=  FACN*1.0 

GO  TO  S 

99  COnIINDF. 
tNU 

3200  F OPT  HAN  DIAGNOSTIC  RESULTS  - F OH  EXPAOT 


NO  ERRORS 
LOAD, Sb 
RUN 


I 
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3 POO  FORTRAN  (P.l)  31/03/iSF> 


SUHhO  UTINE  OAUS  S » IN.  IF  SINOLM*  IPOwFP,  ntT) 

COMMON  SA<13>,  At?H.  TSI.  X17B> 

HEAL  M 

INTF6FH  H.  HP|  . HMAX 
ISI'-N  z I 
MM 1 z N- I 
NP)  z H»l 

IF  sInglh  * n 

TEN  PowEH  * 1 0*  ()••  1-lPOixtKI 

HO  I Pz),  NM) 

AMAX  = AHSF<A(B,P)) 

UMAX  z H 
HP ] = H*I 
00  1 1 zHP  1 * N 
AHSFA  = ABSF  IMIiR)  ) 

IF  I AHSFA  ,(»T.  AMAX)  o 3 
? AM Ax  z AHSFA 
PM 6 X z I 

3 CONTINUE. 

1 F ( AMAX  .LT.  TFN  POhEH)  S.  a 
S IF  SlNGLH  = 1 
A IF  (H  .FU.  nt'AX)  10*  H 
R III)  * l zK  * N 
A T F -IP  s A ( H . I ) 

A ( P . I ) z A (PMA  X « | ) 

A A ( p'  A X . I ) z A TF  "P 
Hit  P = A (H.MP1 ) 

A(H*NP|I  z A(WM.tX*  MP|  | 

A <H,  AX,NP|  | z «(FmP 
IST-N  z - I S I ON 
1,1  u0  1 IzHPl,  N 

M z -,V<1.P)/A(P«P| 

00  1 J*HP|,  N 

7 A(l.J)  z A((.J).  MOA (P* J) 

1 A(I.NPI)  = All*  MP|>  ♦ M*A ( P * NP|> 

IF  (AHSF  (A(N.N)  ) .It.  TENPO»F.P(  ]|,  H 
11  IF  SlNfiLP  = I 
M fONllNUE 

00  21  1*1*  N V 

1 Ml  z 1 _ 1 

NN  * N-IMI 

XX  z A ( NN « NHM 

00  22  J*  1 * I M | 

NP 1 1 z NPI  - J 

22  XX  = XX-A(NN*NPI J)*X(NP1J) 

21  XINN)  z XX/AINN.NN) 

OtT  * \.n 
00  31  I z 1 , n 
31  OET  « OE  T • A ( 1*1) 

I'ET  z 01  T*  I SION 
EN  • 

3200  KIHTPAN  O I AGNOSTIC  PE  SUL  I S - F OH  GAUSSP 


NO  EHHOPS 


59 


i?oo  fiHjioN  <?.n  n/,n/6ft 

SUHwOllTIrF  fNDFXI  JJ*1.J,N,>) 
jj=  jmii-1) • iN?-n  ) /? 

fcNn 

3?00  FOWTWAN  DIAGNOSTIC  WF  SUL  T S - Fmu  InO-X 


NO  EMROWS 
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32U0  FOMTHaN  (2.1)  31/03/66 

SUBROUTINE  HlGSFT  (N) 

COMMON  SA(13),  SB ( 1 2 ) « CI78.  7«»> . X(7a>.  T(12.  12),  E(l?.  12). 

1 A ( 1 2 • 12) • H ( 12«  12) 

NN= (N» IN* 1 ) > /2  * NNP 1 = NN*1 

DO.  7 1*1  . N 
00  7 J*l.  N 

7 H (I  . J)  * E ( I <NI  *E  ( UN) 
h(n<N)=  MN.NI-I.O 
00  I 1*1.  NN 
no  t J*i,  nnp) 

1 C ( 1 « J)  = 0.0 
M?=  2«N  T.  L = o 
00  2 J=l,  N 
JB1=  J*1 

DO  2 1*J.  N 
( *L  » I 
1P1*  1*1 
DO  3 6=1,  ) 

CALL  IN07XIJJ.K.  I.N?) 

< C (L  • JJ>  = C <L . JJ)  . M I .K) 

UO  <i  l\  = JP  I . N 

CALL  iNUF X < JJ.J.A ,N?> 

«.  C (L. JJ) *C (L. JJ) * M l «*) 

1)1)  ->  A*  I . ( 

CALL  I NL)h  X ( J J tK.I.N?) 

A C (L  , JJ)  * C(L.JJ)  ♦ A ( J»K ) 

UO  6 As  IP  1 « N 

CALI.  INUEXIJJ.I.K.M2) 

6 C IL« JJ) =C (L. JJ) *A ( J,M 

2 C(L.NNP| )=  HI  I,  I) 

tNu 

3200  FORTRAN  0 1 AGNUS  TIC  RESULTS  - FOB  HIGSET 


MO  F PROPS 


i?nn  fdwiftan 


(2.1) 


31/01/NN 


SUHWOUI  INF  SOLI  IT  MNIN.M,  I RAW  l TY  I . IP  AWT  T VO  . 1 ST  OW  . I F NWW  1 Ml  , IFXPW|nT) 
COM  .ON  SA(M>«  S*>  ( 1 2 ) « l ( IP)  ,/P  ( 1?)  . *R(1P).  P|N()0>.  X ( SN04  ) , 
t Iw(|P.  IP).  MU.  12).  T(U.  1?) 

IWAwI t r I = » 

1WAW  I t YOsfl 

WIIFFFW  IN  t|.|l(^AXI».  MAX.) 

3ll  1.0  10  Oil.  31.  PP.  3).  UN  I T S TF  ( l I 
31  IwMAX.  M A A * 

HUFFtW  IN ( I . I t ) « I | ) .« ( UMAX) ) 
t 1*0  TO  I).  ?.  P?.  3).  ONITSTFtn 
3 I M Ail  I I Y 1 * I 
11 F 3 . IHN 
? F,M1*  «.,) 
mPPs 
I 1 *N 

I ST  AM  I * •> 

1)0  . 1*1,  N 

/IT)*  ..U 
110  h J*  I , I 

a /in*  r ( n ♦ r ( i . j>  •»  ( j) 

X ( 1 I * il.il 

00  |*|.  MW  1 
J=  iM-'-I 

S X ( I ) * 4 ( ] > .SH ( T ) • x I J) 

7 DO  M I * I ST  AWT . ISTOP 

1 PRINT  = I- 1 

•i  11*  1 1 »N 

IF  (Umax  .lt.  in  1".  11 
III  IFIIl-N  .LQ.  I W‘  i A X 1 u.  13 
1?  JJ=  0 
fiO  TO  I* 

I 3 JJ*  I I-1-.MAX 
KL  = 1 1-N 

00  lb  KJ*1.  JJ 
*L  = KL ♦ I 
IS  a ( K . I)  s w ( K L ) 

14  HUF  F tW  IN  ( I « I ) (MAX  X.  MAX  W) 

3P  1.0  TO  (3W.  13,  34.  3)  « UNlISTFII) 

14  IF(IPwlM)  IS, 

15  IF  (IF  N PHUT  .FT.  1)  TO,  T1 
Til  wwiol  T2,  IWWINT 

I ? FORMAT  (////.  I 10.  30H  POINTS  awf:  I'M  THIS  UaT«  h|  OCk  . /) 

71  IF  (IF  X PRINT  .F  0.  1 ) 73.  74 

73  PRINT  PH,  (X(J).  |=|,  IRRlwl) 

74  bUFFtR  OUT  (2.1)  ( IPR IN  I , IPRImT) 

NP  GO  IT)  (6P.  S3.  S3.  PS),  UnIISTF(P) 
t»  3 HOFFtw  OUI  (P»  1)  ( X ( 1 1 , X(IPrtlNTI) 

H4  GO  10  (t>4.  ?P.  PP.  ?b)  . IINUSTFIPI 

II  UM\x*  MAX  x ♦ J J 
JJ* JJ. I 

hUF  F SR  In  ( 1 , I I (4(  JJ)  , w ( Ixm/iX  ) 1 
XI  1,0  T 0 (P 3.  24.  2P » 3>.  ONITSTF(I) 

24  1I  = ‘ 

1,0  TO  ■ 

11  JJ*  II-N 

1)0  I f>  J*l,  M 
JJ* JJ. 1 

IN  *R(  M*  4 ( J J I 
00  17  K*1 , N 
2/P*0.l> 


00  18  J=l*  N 

|H  U\>x  /ZP.F  <K,.J)  *7  ( J)  ♦ 1 H(K  , J)  ( J) 
17  7P(M=  ZZF> 

XXxtl.l) 

DO  IV  l\  = l.  MP1 
.1=  MP2-K, 

IV  xX=XX.St)(M»ZP(J> 
x ( n = xx 


00  8 J= 1 , N 
• / ( J>  = ZP(  J) 


It  (IF 

N PHIOT  .HO.  1) 

7S . 76 

7S 

PH  I NT 

72,  [STOP 

70 

IF  ( IF 

X PH  1 0 T .tO.  11 

77,  78 

77 

PRToT 

'H.  ( X ( l)  , 1 = 1 , 

1STOP) 

28 

F OpHft [ 

(bE2'l.  101 

78 

HUFFtH 

DOT  17. ! ) ( (STOP. 

ISTOP) 

80 

00  TO 

( bn  * 61.  b 1 , 2b) 

. UNlTSTF (2) 

*1 

HUFFFH 

OIIT  (?» 1 > ( X < 1 > • 

X ( 1 STOP) ) 

21 

00  TO 

(21.  6 . 6 . ^S ) * 

UNlTSTF (2) 

26 

IPAhI  1 Yi)=l 

HF  7 OHO 
6 IST«HI=| 

00  10  1 
2?  1-NO  FILE  2 

continue 

l-.N' ' 


Von  F OH  THAN  III*'  HOST  ic  HFSOlTS  - F''P  SOLOMON 


UNDEFINED  SIMPLE.  MXHHI.1S 

MAXl/ 

FQUIP.  2*h  T 

LOAD.Sb 

PUN 


APPENDIX 


I 

| I 

l 

In  addition  to  the  "white  sequence,"  { e 1 ) , 1 5 i § 4266,  discussed 
in  the  text  of  this  paper,  seven  other  "white  sequences"  have  been  generated 
I by  IELASHO.  These  eight  "white  sequences"  have  been  permanently  stored  on 

| magnetic  tape;  this  tape  is  labeled  Tape  B. 

The  eight  "white  sequences"  (in  the  order  in  which  they  are  stored 
| on  Tape  B)  are: 

v[L)  = (01),  ins  4oo4;  e = 9.63778^51  * i°g10N 

K„  ■ v[2)  = {e1},  1 s i s 4oio;  e = 1.27323954  * V" 

I vp-1  = (91),  1 § i * U006;  0 = 3.0U800610  * 0.2( number  of  cm/ft) 

| V[M  = (01),  1 5 i § kO06 ; 9 = l.U6^59l89 

P]  - (01),  1 S i S 4006;  0 = 4.9711*9872  * log10it 
■ v[6)  = {01},  1 § i § 4266;  0 = 2.71828183  * e 

I v<7)  = (01),  1 5 i 5 4008;  e = 2.30258509  * 1/N 

v<8>  = (e1) , 1 § i s U006;  0 = 5.72957795  * o.i(  180/n) . 


f- 


<1 


l 


6^ 


I DEC 

t(8) 


(8) 


IX<8>(L<8)> 


IX(8)(1)  • • 

end-of-file 

end-of-file. 

Here  IDEC,  L,  and  IX(l)  are  defined  the  same  as  in  the  section  of 
the  text  on  ISLASHO. 

The  "-white  sequence"  on  Tape  B can  he  used  as  direct  input 


,(s) 


to  GAUSSIAN.  In  order  to  use  , 2 § i % 8,  as  input  to  GAUSSIAN,  either 


the  specified  "white  sequence"  must  first  he  placed  on  an  auxiliary  tape  or 
GAUSSIAN  must  he  modified  to  skip  over  previous  "white  sequences." 

The  first  UOOO  points  of  each  of  the  eight  above  "white  sequences" 


have  heen  interlaced  hy  the  scheme 

(8) 

' * V8( i-l)  + 8 ’ vi  ' 
for  1 § i 5 UOOO.  This  32,000  point  "white  sequence"  has  also  heen  stored 


v - v(l)  v - v(2) 

V8( i-l)  + 1 - Vi  ' V8( i-l)  + 2 - Vi  » 


permanently  on  a second  magnetic  tape  which  is  labeled  Tape  C. 
Tape  C has  the  format : 

500 


500 


500 
v. 


31,501  32,000 

end-of-file 
end-of-file. 

Tape  C can  be  used  directly  as  input  to  GAUSSIAN. 


i :r 


'I WgKlfy 
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Tape  B has  the  following  format: 


500 

. 

3501 

k 


v!1} 

Uooo 


end-of-file 


1X^(1)' 


end-of-file 


(8)  (8) 

V^00l"  * Vl+00* 

end-of-file 


REFERENCES 


J.  N.  Franklin,  "Numerical  Simulation  of  Stationary  and  Non- 
stationary Gaussian  Random  Processes,"  SIAM  Review,  No.  1,  68 
( J anuary,  1965 ) . 

Ibid. , pp.  68-80. 

J.  S.  Bendat,  Principles  and  Applications  of  Random  Noise  Theory 
(John  Wiley  and  Sons7"NewYork7^^58Tr-pF^^OiT^! 

k Ibid.,  pp.  lUl,  161-3. 

<5 

R.  B.  Blackman  and  J.  W.  Tukey,  The  Measurement  of  Power  Spectra 
(Dover  Publications,  Inc.,  New  York,  195^),  PP-  35-7- 

^ Franklin,  0£.  cit . , pp.  68-80. 

7 

W.  B.  Davenport  and  W.  L.  Root,  An  Introduction  to  the  Theory  of 
Random  Signals  and  Noise  (McGraw-Hill,  New  York,  1958),  pp.  182-3, 
225-7. 


8 


Ibid.,  pp.  232-4,  375-6. 


J.  N.  Franklin,  "Deterministic  Simulation  of  Random  Processes," 
Math.  Comp.,  _1£,  47-50  (1963). 

^ G.  E.  P.  Box  and  Mervin  E.  Muller,  "A  Note  on  the  Generation  of 
Random  Normal  Deviates,"  Ann.  Math.  Statist.,  29,  610-1  (1958). 

^ J.  N.  Franklin,  "The  Covariance  Matrix  of  a Continuous  Auto- 
regressive Vector  Time  Series,"  Ann.  Math.  Statist.,  34,  1259-64 
(1963). 


12 

F.  B.  Hildebrand,  Introduction  to  Numerical  Analysis  (McGraw-Hill, 
New  York,  1956),  pp.  486-9. 

^ L.  Fox,  An  Introduction  to  Numerical  Linear  Algebra  (Oxford 
University  Press,  New  York,  1965),  PP-  60-91. 

14 

J.  N.  Franklin,  "Numerical  Simulation  of  Stationary  and 
Nonstationary  Gaussian  Random  Processes,"  SIAM  Review,  I, 

No.  1,  70  (January,  1965). 

15 

Defense  Research  Laboratory,  The  University  of  Texas,  Acoustical 
Report  No.  252  resulting  from  research  under  U.  S.  Naval  Oceanographic 
Office  Contract  N62306-I358,  8 December  1965. 

66 


Defense  Research  Laboratory,  The  University 
Report  No.  1 on  Contract  DA  28-043  AMC-003l6(E), 
1965,  pp.  27,  31. 


• > 


of  Texas,  Annual 
1 July  1964-30  June 


17 


Ibid 


P.  23. 


BIBLIOGRAPHY 


Bendat,  J.  S.,  Principles  and  Applications  of  Random  Noise  Theory 
(John  Wiley  and  Sons,  New  York,  1956 ) . 

Blackman,  R.  B.  and  J.  W.  Tukey,  The  Measurement  of  Power  Spectra 
(Dover  Publications,  Inc.,  New  York,  1956) • 

Box,  G.  E.  P.  and  Mervin  E.  Muller,  "A  Note  on  the  Generation  of 
Random  Normal  Deviates,”  Ann.  Math.  Statist.  29 > 6l0-l(l958). 

Davenport,  W.  B.  and  W.  L.  Root,  nn  Introduction  to  the  Theory  of 
Random  Signals  and  Noise  (McGraw-Hill,  New  York,  195®) • 

Ellis,  G.  E.  and  R.  L.  Boston,  "Program  SPECT:  Power  Spectral 

Analysis  of  a Random  Process  Using  the  CDC  1604  Computer," 
Defense  Research  Laboratory,  The  University  of  Texas, 

Austin,  Texas,  Acoustical  Report  A-252  (December,  1965). 

Fox,  L.,  An  Introduction  to  Numerical  Linear  Algebra  (Oxford 
University  Press,  New  York,  19$?) • 

Franklin,  J.  N.,  "The  Covariance  Matrix  of  a Continuous  Auto- 
regressive Vector  Time  Series,"  Ann.  Math.  Statist, 
lit,  1259-64  (1965). 

Franklin,  J.  N.,  "Deterministic  Simulation  of  Random  Processes," 
Math.  Comp.  1J,  28-59  (1965). 

Franklin,  J.  N.,  "Numerical  Simulation  of  Stationary  and  Non- 
s ationary  Gaussian  Random  Processes,"  SIAM  Review  1, 

No.  1,  68-00  (January,  196?). 

Hildebrand,  F.  B.,  Introduction  to  Numerical  Analysis  (McGraw-Hill, 
New  York,  195$7» 

Richie,  W.  C.  and  D.  R.  Chick,  "Report  on  Research  to  Investigate 
Atmospheric  Infrasound,"  Defense  Research  Laboratory, 

The  University  of  Texas,  Austin,  Texas,  Acoustical  Report  A-249 
(November,  1965). 


68 


MR6MMWB  I ■ 


r 


27  October  1966 


DISTRIBUTION  LIST  FOR  DRL-A-258 
UNDER  CONTRACT  NObsr- 93175 


Copy  No. 

1 - 4-  Commander,  Naval  Ship  Systems  Command 

Department  of  the  Navy 

Washington,  D.  C.  20360 

Attn:  Mr.  C.  D.  Smith  (Code  1622C ) 

5 Commanding  Officer 

U.  S.  Navy  Mine  Defense  Laboratory 
Panama  City,  Florida  32402 
Attn:  Mr.  Carl  M.  Bennett 

6 Commanding  Officer 

U.  S.  Naval  Oceanographic  Office 

Suitland,  Maryland  20390 

Attn:  Mr.  R.  C.  Guthrie  (Code  7240) 

Oceanographic  Analysis  Division 
Marine  Sciences  Department 

7 Commanding  Officer 

U.  S.  Navy  Electronics  Laboratory 
San  Diego,  California  92152 
Attn:  Code  3180 

8 Office  of  Naval  Research 
Resident  Representative 
The  University  of  Texas 
2507  Main  Building 
Austin,  Texas  78712 

9 Acoustics  Division,  DRL/UT 


10 

Signal  Physics  Branch, 

11 

K.  J. 

Diercks,  DRL/UT 

12 

G.  E. 

Ellis,  DRL/UT 

15 

D.  W. 

Everts on,  DRL/UT 

14 

K.  W. 

Harvel,  DRL/UT 

15 

S.  P. 

Pitt,  DRL/UT 

16  - 17 

Library,  DRL/UT 

